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ICE  DYNAMICS 


William  D.  Hibler  III 


INTRODUCTION 

Sea  Ice  dynamics  deals  with  the  way  momentum  is  transferred  through 
the  sea/lce  system.  The  essential  features  of  this  transfer  may  be  charac¬ 
terized  by  equations  representing  the  following  elements:  1)  a  momentum 
balance  describing  ice  drift,  including  air  and  water  stresses,  Coriolis 
force,  internal  ice  stress,  inertial  forces,  and  ocean  current  effects;  2) 
an  ice  rheology  which  relates  the  ice  stress  to  the  ice  deformation  and 
strength;  and  3)  an  ice  strength  determined  primarily  as  a  function  of  the 
ice  thickness  distribution.  In  discussing  these  components  it  is  implicit¬ 
ly  assumed  that  a  description  of  sea  ice  as  a  two-dimensional  continuum  is 
possible.  Such  a  continuum  description  does  not,  however,  rule  out  using 
these  equations  to  represent  rapidly  changing  ice  conditions  as  occur,  for 
example,  near  the  ice  margin.  This  can  be  accomplished  by  allowing  various 
quantities  (e.g.  ice  mass,  ice  strength,  water  stress)  to  vary  rapidly  in 
space.  Also,  most  present-day  treatments  of  ice  dynamics  employ  highly 
nonlinear  rheologies  which  can,  as  will  be  demonstrated  later,  introduce 
discontinuities  into  the  ice  velocity  field,  even  though  the  forcing  and 
mass  fields  may  be  rather  smooth. 

In  the  following  review,  these  three  elements  are  discussed  in  order. 
In  this  discussion  the  emphasis  is  on  understanding  the  effect  of  ice  rhe¬ 
ology  on  ice  drift  and  deformation,  and  both  interior  and  ice  edge  dynamics 
are  examined.  The  final  section  examines  the  role  of  ice  dynamics  in  ice- 
atmosphere-ocean  coupling. 


MOMENTUM  BALANCE 

The  momentum  balance  describes  the  forces  that  determine  drift  and  de¬ 
formation  of  sea  ice.  Denoting  the  ice  mass  per  unit  area  by  m,  the  momen¬ 
tum  balance  in  a  two-dimensional  Cartesian  coordinate  system  is 

mDtu»-mfkxu+f  +  t  +  F-  mgVH  ( 1 ) 

t  —  —  —  —a  — w  — 

where  D^  ■  (3/3t  +  •  V)  is  the  substantial  time  derivative,  k  is  a  unit 

vector  normal  to  the  surface,  ii  is  the  ice  velocity,  f  is  the  Coriolis  pa¬ 
rameter,  _ra  and  are  the  forces  due  to  air  and  water  stresses,  H  is  the 
the  dynamic  height  of  the  sea  surface,  and  £  is  the  force  due  to  variations 
in  internal  ice  stress.  In  this  formulation,  Jj,  is  assumed  to  include 
frictional  drag  due  to  the  relative  movement  between  the  ice  and  the  under- 


lying  ocean.  Also,  H  can,  in  principle,  include  the  variation  of  the  sea 
surface  height  due  to  atmospheric  pressure  changes  as  well  as  to  geostroph- 
ic  current  balance.  In  this  momentum  balance,  the  air  and  water  stresses 
are  normally  determined  from  idealized  boundary  layers,  assuming  constant 
turning  angles  (McPhee  1979,  Brown  1980,  Leavitt  1980): 

t  =  c'  (U  cos<t>  +  k  x  U  sin4>)  (2) 

—a  a  -g  —  — g 

t  -  c'[(U  -  u)  cos0  +  k  x  (U  -  u)  sin9)]  (3) 

— w  w  — w  —  —  — w  — 


where  Ug  is  the  geostrophic  wind,  Uj,  the  geostrophic  ocean  current,  ca  and 
c^  air  and  water  drag  constants,  and  <f>  and  6  air  and  water  turning  angles. 
In  practice  both  the  currents  and  the  sea  surface  tilt  are  estimated  from 
geostrophic  considerations  by  setting  H  equal  to  the  dynamic  heights  and 
computing  currents  by  Uw  *  gf-1k  x  7h,  where  g  is  the  acceleration  due  to 
gravity.  In  general,  both  ca  and  c^  are  nonlinear  functions  of  the  winds 
winds  and  currents.  The  two  most  commonly  used  formulations  are  1)  linear, 
where  ca  and  c^  are  taken  to  be  constant,  and  2)  quadratic,  where 


c' 

a 


(4) 


and 


c' 

w 


P  c 
w  w 


(5) 


with  ca  and  being  dimensionless  drag  coefficients  (with  typical  values 
of  0.0012  and  0.0055  respectively;  McPhee  1980). 


There  are  also  a  variety  of  other  nonlinear  boundary  layer  formula¬ 
tions,  several  of  which  have  been  reviewed  by  McPhee  (1982).  In  some  of 
these  boundary  layer  models,  a  precise  drag  law  cannot  be  specified  in  the 
form  of  eq  3  and  5,  but  rather  the  boundary  layer  characteristics  are  de¬ 
termined  by  the  stress  at  the  ice/ water  interface. 


Basically,  in  this  formulation  the  essential  causes  of  the  ice  motion 
can  be  thought  of  as  the  geostrophic  wind  above  the  atmospheric  boundary 
layer  and  the  ocean  current  beneath  the  oceanic  boundary  layer.  These 
forces  are  transmitted  to  the  ice  via  simple  integral  boundary  layers.  In 
addition,  a  steady  ocean  current  introduces  a  tilt  (7h)  of  the  sea  surface 
height  (or  alternatively  a  pressure  on  the  ice),  which  also  affects  the  Ice 
motion.  A  subtle  but  important  point  with  regard  to  the  oceanic  boundary 
layer  is  that  there  is  no  stress  transferred  from  the  bottom  of  the  bounda¬ 
ry  layer  to  the  deeper  ocean.  Instead,  the  drag  upon  the  ice  arises  be¬ 
cause  of  a  combination  of  shearing  and  rotation  in  the  boundary  layer. 

This  point  has  important  ramifications  for  ice-ocean  coupling.  Specifical¬ 
ly,  while  the  stress  transferred  to  the  ocean  is  equal  to  the  wind  stress 
less  the  ice  interaction,  it  will  in  general  be  somewhat  different  than  the 
water  drag  on  the  ice.  This  can  be  seen  by  taking  the  curl  of  eq  1  with 
the  Internal  ice  force  tqrra  set  equal  to  zero.  When  this  is  done  the  curl 
of  the  water  drag  will  differ  from  the  curl  of  the  wind  stress  by  a  term 
proportional  to  mfVm  which  arises  because  of  the  Coriolis  force  on  the 
ice. 
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Relative  magnitude  of  forces 


With  respect  to  the  dominant  components  of  the  momentum  balance,  ob¬ 
servations  in  conjunction  with  dimensional  analysis  show  that  the  force 
balance  is  largely  dominated  by  the  Coriolis  term,  the  air  and  water 
stresses,  and  the  ice  interaction.  Specifically,  measurements  of  water  and 
air  stresses  show  them  to  be  of  the  order  of  0.1  N  m“  .  For  3-m-thick  ice 
moving  at  0.1  m  s'*  ,  the  Coriolis  force  is  about  0.05  N  m~  ,  a  value  com¬ 
mensurate  with  the  wind  and  water  stresses.  With  respect  to  inertial 
force,  for  3-m-tljick  ice,  the  ice  must  accelerate  up  to  a  typical  drift 
rate^f  0.1  m  s“  in  less  than  15  minutes  to  yield  an  inertial  force  of  0.1 
N  m“  .  However,  the  acceleration  term  can  be  significant  in  Inertial  os¬ 
cillations.  In  this  case  an  acceleration  is  induced  by  changing  the  direc¬ 
tion  of  the  ice  motion  on  the  same  time  scale  as  the  inverse  of  the  Corio¬ 
lis  parameter  (about  3  hours).  For  most  applications  the  momentum  advec- 
tion  term  is  very  small,  since ^elocity  changes  of  0.1  m  s~  per  kilometer 
would  be  needed  for  a  0.1  N  m~  force. 

With  respect  to  geostrophic  current  effects  and  the  ocean  tilt  terms, 
their  importance  depends  on  the  magnitude  of  the  current  fields,  which  are 
generally  smaller  than  the  ice  drift  rates.  In  the  central  Arctic,  for 
example,  empirical  model  studies  (Hibler  and  Tucker  1979)  indicate  that 
steady  current  effects  are  normally  quite  small  (several  percent  relative 
to  the  wind-driven  components)  for  ice  velocities  averaged  over  a  few  days. 
However,  even  in  the  central  basin,  current  and  tilt  effects  are  signifi¬ 
cant  over  a  long  time  period  since  the  wind  effects  tend  to  average  out 
(see  Fig.  1).  Also,  in  relatively  shallow  regions,  where  barotropic  ef¬ 
fects  are  significant,  the  role  of  currents  in  ice  drift  may  be  critical. 

The  remaining  term  in  the  momentum  balance  is  the  ice  interaction. 
Under  typical  winter  conditions  in  the  central  Arctic  this  can  be  shown  to 
be  comparable  to  the  other  terms  by  residual  calculations.  A  force  balance 
obtained  by  Hunkins  (1975)  employing  such  a  procedure  is  shown  in  Figure  2. 
A  somewhat  more  direct  method  is  to  make  use  of  the  fact  that  in  enclosed 
areas,  such  as  the  Bay  of  Bothnia,  the  ice  does  not  move  at  all,  even  with 
significant  winds  (Lepparanta  1980,  Omstedt  and  Sahlberg  1977).  Based  on 
the  wind  drag  values  discussed  above  this  indicates  ice  stresses  at  least 
of  the  same  order  of  magnitude  as  the  wind  and  water  drag.  It  is  also  pos¬ 
sible  to  estimate  ice  stresses  by  considering  the  mechanics  of  the  ridging 
process.  More  details  on  this  approach  and  the  relation  of  the  ice  inter¬ 
action  to  the  ridging  process  will  be  discussed  in  the  section  beginning  on 
p.  23. 

Free  drift  analysis 

Much  of  this  paper  will  be  concerned  with  analyzing  the  effects  and 
characteristics  of  the  internal  ice  stress  term  F.  However,  it  is  instruc¬ 
tive  to  first  examine  a  special  case  of  the  momentum  balance  in  which  there 
is  no  ice  Interaction.  This  case  is  usually  referred  to  as  free  drift.  To 
further  simplify  matters,  only  steady  free  drift  will  be  considered,  so 
that  the  inertial  terms  can  also  be  removed. 

Following  McPhee  (1980)  and  Thorndike  and  Colony  (1982)  it  is  conve¬ 
nient  to  use  complex  notation  to  analyze  this  case.  Writing  all  vectors  in 
complex  form,  i.e.  A  -  AeiY,  and  noting  that  the  force  on  the  ice  due  to 


Figure  1.  Simulated  and  observed  net  drift  of  three  ice  stations 
over  a  two-year  period  employing  an  eopirical  linear  ice  drift 
model.  Simulations  are  shown  both  with  and  without  current  effects 
(from  Hibler  and  Tucker  1979). 
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Figure  2.  An  estimate  of  the  force 
balance  on  sea  ice  for  winter  condi¬ 
tions  based  on  wind  and  water  stress 
measurements  (from  Hunkins  1975).  In 
this  balance  the  force  due  to  inter¬ 
nal  ice  stress  is  determined  as  a  re¬ 
sidual  and  the  dashed  line  shows  the 
ice  velocity. 


sea  surface  tilt  is  given  by  -imfjj,,,  the  steady  state  momentum  balance  can 
be  written  in  the  form 


i  <♦>  19 

p  c  U  U  e  *  «  imfU  +  p  c  U  U  e 
a  a~g  g  ~  w  w~  — 


(6) 


where  U  is  the  ice  velocity  relative  to  the  steady  currents: 


U  «  U  -  U. 


(7) 


Letting  <5  be  the  clockwise  angle  from  geostrophlc  wind  to  ice  drift  direc¬ 
tion,  eq  6  can  be  rewritten 


„2  i«  „2  -i<f> 
U  e  =  U  e 
g 


P  c 
w  w 

P  c 
a  a 


i  9 


(e~ '  + 


mf 


iff/ 2 


P  c  U 
w  w 


)• 


(8) 


Several  things  are  clear  by  inspec¬ 
tion  of  eq  8.  For  one,  if  the  ice 
mass  is  small  enough  to  make  the 
Coriolis  term  negligible,  then  the 
ice  will  drift  at  a  rate  equal  to  a 
fixed  fraction  of  the  geostrophic 
wind  with  a  drift  angle  0-$  rela¬ 
tive  to  the  geostrophic  wind.  As 
the  mass  increases  the  ice  will 
tend  to  drift  further  to  the  right 
of  the  geostrophic  wind.  Similar 
arguments  also  apply  to  the  ice 
speed,  with  larger  ice  speeds  caus¬ 
ing  the  ice  to  drift  more  nearly 
parallel  to  the  geostrophic  wind. 

A  final  point  of  note  is  (referring 
back  to  eq  6)  that  for  linear  boun¬ 
dary  layers,  the  ratio  and  angle 
between  the  ice  drift  and  the  wind 
will  be  fixed  for  a  given  ice  mass, 
with  increasing  ice  mass  causing 
the  ice  to  drift  further  to  the 
right  of  the  wind. 

A  more  quantitative  illustra¬ 
tion  of  the  nonlinear  boundary  lay¬ 
er  case  is  given  in  Figure  3  which 
shows  solutions  of  eq  8  in  the  form 
of  ice  speed  and  turning  angle  ver¬ 
sus  geostrophic  wind  speed  for  dif¬ 
ferent  ice  masses.  In  this  figure 
drag  coefficients  and  Coriolis  val¬ 
ues  similar  to  those  used  by  Hibler 
(1979)  were  employed  (pw  ■  10 3  kg 
m-  *  Pa  *  1.3  kg  m~  ,  Cy  *  0.0055, 


Figure  3.  The  ice  drift  rate  versus 
the  geostrophic  wind  speed,  and  the 
clockwise  turning  angle  from  the  geo¬ 
strophic  wind  to  the  ice  velocity. 
These  curves  were  calculated  assuming 
free  drift,  and  quadratic  boundary 
layers.  Curves  for  different  ice 
thicknesses  are  shown. 


ca  -  0.0012,  f  -  1.46  x  10"  s~  ,  6  ■  $  -  25°).  One  interesting  result 
of  this  plot  is  the  very  nearly  linear  relation  between  the  wind  speed  and 
ice  drift  for  all  but  very  small  speeds,  despite  the  fact  that  boundary 
layer  drags  are  nonlinear.  This  situation  arises  because  both  boundary 
layers  are  quadratic,  and  hence,  except  at  low  speeds,  largely  balance  one 
another.  However,  the  turning  angle  does  show  somewhat  greater  sensitivity 
to  wind  speed.  Also  note  that  although  U  and  Ug  are  linearly  related, 
the  best  linear  fit  does  not  go  through  the  origin. 

Effects  of  ice  interaction;  A  qualitative  discussion 

It  is  possible  to  get  a  qualitative  idea  of  the  effect  of  ice  inter¬ 
action  by  examining  the  free  drift  force  balance.  The  most  obvious  effect 
of  the  ice  stress  is  to  change  the  ice  drift  direction  and  magnitude  from 
the  free  drift  case.  While  there  can  be  exceptions,  in  many  cases  the  net 
effect  of  the  ice  interaction  on  the  local  force  balance  is  to  cause  a 
force  opposing  the  wind  stress  and  the  Coriolis  force,  roughly  in  the  man¬ 
ner  shown  in  Figure  4.  As  a  consequence,  to  achieve  the  same  ice  velocity 
under  these  circumstances  a  larger  wind  stress  more  nearly  parallel  to  the 
ice  velocity  is  required. 

Less  obvious  effects  arise  from  the  nonlinear  nature  of  the  ice  inter¬ 
action.  As  will  be  discussed  later,  over  a  wide  range  of  strain  rates  the 
ice  stress  is  expected  to  change  relatively  little.  As  a  consequence  the 
force  due  to  ice  stress  is  less  sensitive  to  ice  speed.  Because  of  this, 
the  modifications  to  free  drift  will  become  particularly  pronounced  at  low 
wind  speeds  and  ice  velocities.  Nonlinear  interaction  can  also  give  rise 
to  the  interesting  situation  where,  although  the  ice  is  strongly  interact¬ 
ing,  its  local  velocity  may  obey  free  drift.  This  arises  because  the  force 
due  to  internal  ice  stress  is  given  by  the  gradient  of  the  ice  stress,  so 
that  even  though  the  stress  is  high  its  gradient  may  be  very  small.  Be¬ 
cause  of  complexities  like  this,  high  correlation  of  ice  motion  with  wind 
fields  cannot  be,  by  itself,  taken  as  evidence  for  the  absence  of  ice  in¬ 
teraction. 

A  third  subtle  situation  which 
can  arise  due  to  nonlinear  ice  in 
teraction  is  a  much  more  "rigid 
body"  motion  of  large  regions  of 
sea  ice.  Within  the  framework  of 
a  plastic  rheology  this  state  can 
arise  where  the  ice  stress  is  not 
large  enough  to  exceed  the  plastic 
yield  value.  Under  this  situation 
the  motion  of  the  rigid  body  will 
behave  very  similarly  to  a  spatial-  rw  mtu 

ly  "averaged"  free  drift,  with  the 

main  difference  being  that  the  de-  Figure  4.  If  the  wind  stress  Ta 

formation  will  be  much  different  balances  the  water  drag,  Coriolis 

than  what  would  be  expected  from  and  current  forces,  then  a  larger 

free  drift  considerations.  This  stress  t^,  turned  further  to  the 

particular  situation  is  discussed  right,  is  required  to  balance  the 

in  more  detail  in  a  later  section.  forces  with  the  addition  of  a  force 

F  due  to  internal  ice  stress. 


As  these  general  remarks  Indicate,  the  effects  of  ice  interaction  can 
manifest  themselves  in  subtle  and  sometimes  counter-intuitive  ways.  In  the 
remainder  of  this  chapter  some  of  these  subtleties  will  be  examined  in  more 
detail. 


SEA  ICE  RHEOLOGY 

Analysis  of  the  momentum  balance  has  shown  the  importance  of  the  in¬ 
ternal  ice  stress  in  ice  drift.  Because  of  this,  the  determination  of  a 
realistic  constitutive  law  relating  the  ice  stress  to  the  deformation  is 
critical  to  understanding  ice  drift  and  deformation.  Some  of  the  main 
questions  relating  to  ice  rheologies  are  a)  the  degree  of  nonlinear  versus 
linear  behavior  for  both  shearing  and  compressive  failure,  b)  the  relative 
magnitudes  of  shear  and  compressive  strengths,  c)  the  presence  (if  any)  of 
an  unconfined  pressure  term  due  to  random  bumping  of  floes,  and  d)  the  mag¬ 
nitude  and  dependence  of  ice  strength  on  ice  thickness  and  floe  size  char¬ 
acteristics.  These  questions  are  important  for  understanding  the  way  the 
ice  responds  to  external  forcing  and  for  understanding  the  behavior  of  the 
coupled  air-ice-ocean  system.  In  this  section  observational  and  theoreti¬ 
cal  work  relating  to  the  first  three  questions  is  discussed.  The  depen¬ 
dence  and  relationship  of  the  ice  strength  to  thickness  characteristics  is 
discussed  in  the  next  section.  An  attempt  is  made  to  present  first  a  fair¬ 
ly  unified  general  framework,  and  then  to  expand  on  the  framework  with  var¬ 
ious  pedagogical  examples  employing  analytical,  numerical  and  observational 
work. 

In  discussing  the  rheology  of  sea  ice  it  will  be  assumed  that  a  de¬ 
scription  of  sea  ice  as  a  two-dimensional  continuum  is  possible.  While 
such  an  approximation  is  reasonable,  it  should  be  remembered  that  pack  ice 
does  consist  of  individual  floes,  which  can  be  as  large  as  20  km  in  the  in¬ 
terior  of  the  Arctic  Basin  and  only  tens  of  meters  across  near  the  ice  mar¬ 
gins.  Because  of  these  considerations  the  scale  at  which  the  continuum  ap¬ 
proximation  is  valid  will  vary.  More  generally,  continuum  motion  has  to  be 
considered  a  statistically  defined  quantity  for  sea  ice. 

Theoretical  framework 

Considering  sea  ice  to  be  a  two-dimensional  isotropic  continuum,  the 
most  general  constitutive  law  (see,  e.g. ,  Glen  1970,  Malvern  1969)  can  be 
written  in  the  following  form: 

°ij  “  2n^ij  +  -  P]  6^  (9) 

where  ojj  is  the  two-dimensional  stress  tensor  and  the  two-dimensional 
strain  rate,  *  1  for  i  equal  to  j  or  0  for  i  not  equal  to  j  and  n  and  £ 

can  vary  spatially.  In  eq  9,  n,  C  and  P  are  in  general  functions  of  the 
two  invariants  of  the  strain  rate  tensor.  These  invariants  can  be  taken  as 
the  principal  values  of  ice  strain  rate  tensor,  or  alternatively  as  and 
fcjj  ^ij*  l  In  these  expressions  and  equations  (and  throughout  this  Mono¬ 
graph)  repeated  subscripts  are  summed  over.  For  example  6^  ■  +  £yy]« 

Within  the  framework  of  this  equation,  n  and  C  can  be  thought  of  as  non¬ 
linear  shear  and  bulk  viscosities,  while  P  can  be  viewed  as  a  pressure 
term.  The  motivation  for  this  nomenclature  is  the  fact  that  the  compres- 


sive  stresses  and  deviatoric  stresses  can  be  written  as  (as  easily  shown  by 
manipulation  of  eq  9): 
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where  /2  is  the  deviatoric  strain  rate.  To  obtain  an  inter¬ 

nal  ice  force  for  the  momentum  balance,  the  stress  tensor  is  differenti¬ 
ated: 


(12) 


Particular  forms  of  this  constitutive  law  include  both  viscous  and 
plastic  rheologies.  In  general,  a  “viscous"  rheology  is  defined  as  one  in 
which  stress  can  only  be  sustained  through  non-recoverable  dissipation  of 
energy  by  deformation.  A  "plastic"  rheology,  on  the  other  hand,  is  one  in 
which  stress  can  be  sustained  through  lack  of  deformation  or  deformation 
where  the  energy  is  recoverable  (i.e.  elastic  deformation;  e.g.  Malvern 
1969).  Special  cases  of  these  rheologies  often  used  in  sea  ice  are  a 
"linear  viscous"  rheology,  where  the  stress  depends  linearly  on  the  strain 
rates,  and  an  "ideal  rigid  plastic"  rheology,  where  the  stress  state  is 
either  indeterminate  or  independent  of  the  magnitude  of  the  strain  rates. 

In  the  linear  viscous  case,  n,  S  and  P  are  constant.  This  results  in  a 
stress  state  linearly  dependent  on  the  strain  rates.  In  the  rigid  plastic 
case  the  stress  state  is  taken  to  be  fixed,  independent  of  the  magnitude  of 
the  strain  rate  invariants,  provided  the  ratio  of  the  invariants  does  not 
change.  In  this  case  n  and  S  will  be  nonlinear  functions  of  the  invari¬ 
ants.  In  the  rigid  plastic  case,  how¬ 
ever,  ti  and  C  are  not  well  defined  for 
zero  strain  rates.  Under  these  condi¬ 
tions  the  ice  is  understood  to  move 
rigidly,  with  the  stress  determined  by 
external  balances. 


/  Bulk  Viscosity 
Only 


Plastic 


Figure  5.  Allowable  stress 
states  for  a  linear  viscous 
rheology  with  either  a  bulk 
or  shear  viscosity,  and  for  an 
ideal  rigid  plastic  rheology 
with  an  elliptical  yield  curve. 
The  stress  states  are  plotted 
as  a  function  of  the  principal 
components  of  the  two-dimen¬ 
sional  8  tress  tensor. 


Some  of  the  linear  viscous  and 
plastic  rheologies  proposed  for  sea  ice 
are  shown  in  Figure  5,  where  the  allow¬ 
able  stress  states  for  two  particular 
linear  viscous  cases  of  eq  9  are  shown. 
These  states  are  a  viscous  fluid  with 
only  a  shear  viscosity  (n  *  constant,  C 
■  P  =  0)  ,  and  a  fluid  having  resistance 
to  dilation  and  compression  but  not  to 
shear  (C  =  constant,  n  ■  P  *  0) .  For 
an  ideal  plastic  constitutive  law  the 
stress  state  lies  on  or  inside  some 
fixed  yield  curve.  In  Figure  5  an 
elliptical  yield  curve  is  shown.  For  a 
plastic  rheology  a  rule  is  needed  to 
uniquely  relate  the  stress  state  to  the 


strain  rates.  The  most  common  rule  (e.g.  Malvern  1969)  is  to  take  the 
ratio  of  the  strain  rates  for  a  given  stress  state  to  be  that  of  a  vector 
normal  to  the  surface  (this  is  referred  to  as  the  normal  flow  rule).  For 
the  ellipse,  the  normal  flow  rule  yields  n,  £  and  P  values  of  (see,  e.g., 
Hibler  1977) 


P  - 

C  - 

n  - 

where 


P*/2 

(13) 

P*/2  A 

(14) 

S/e2 

(15) 

[( + 

xx 


e  ) 

yy 


[l+(l/e;]  + 


4 

71 


e  +  2e 
xy  xx 


(16) 


-P*  is  the  maximum  compressive  stress  allowable  along  any  principal  stress 
axis,  and  e  is  the  ratio  of  the  lengths  of  principal  axes  of  the  ellipse. 
For  rigid  motion,  the  stress  state  will  lie  within  the  ellipse.  Note  that 
both  linear  viscous  laws  allow  large  amounts  of  tensile  strength  (i.e. 
positive  stresses),  whereas  this  particular  plastic  law  allows  only  small 
tensile  strengths. 


A  general  class  of  rheologies  almost  equivalent  to  the  most  general 
case  (eq  9)  has  been  proposed  by  Hibler  (1979)  under  the  name  “viscous- 
plastic.  "  The  rheology  refers  to  the  case  of  eq  9  in  which  rigid  motion 
does  not  occur,  or,  stated  differently,  any  case  of  eq  9  in  which  knowledge 
of  the  strain  rate  uniquely  specifies  the  stress  state.  This  class  of 
rheologies  was  proposed  partly  for  numerical  reasons  and  partly  based  on 
physical  arguments.  On  the  numerical  side,  the  viscous-plastic  concept  is 
an  efficient  means  of  modeling  highly  nonlinear  sea  ice  dynamics  without 
such  artifacts  as  elastic  waves  which  arise  in  elastic-plastic  simulations 
(Colony  and  Rothrock  1980).  Also,  the  viscous-plastic  approach  lends  it¬ 
self  well  to  formulation  in  a  fixed  Eulerian  grid,  which  is  essential  for 
long-term  seasonal  simulations.  On  the  physical  side  the  viscous  plastic 
rheology  includes  a  more  general  class  of  constitutive  laws  than  can  be 
considered  employing  elastic  plastic  rheologies  (see,  e.g.,  Pritchard 
1975).  Specifically,  with  the  viscous  plastic  approach  it  is  not  necessary 
to  employ  the  associated  flow  rule.  Some  of  these  more  general  cases  may 
be  particularly  applicable  to  the  marginal  sea  ice  zone. 

Linear  versus  nonlinear  behavior 


A  particular  problem  in  early  developments  in  ice  dynamics  was  the 
degree  of  linear  versus  nonlinear  behavior  of  the  ice  rheology.  Early  ef¬ 
forts  to  approximate  the  ice  interaction  made  use  of  linear  viscous  models. 
Laikhtman  (1964),  for  example,  proposed  a  linear  Newtonian  viscous  model 
employing  only  a  shear  viscosity,  which  was  easy  to  deal  with  numerically 
and  was  subsequently  used  by  a  number  of  authors  (Egorov  1970,  Doronin 
1970)  to  carry  out  empirical-numerical  studies.  The  most  well  known  appli¬ 
cation  of  this  rheology  was  a  calculation  by  Campbell  (1965)  of  the  steady 
circulation  of  the  Arctic  Basin  using  an  idealized  mean  annual  wind  field. 
These  empirical  model  studies  demonstrated  that  the  Newtonian  viscous  rhe¬ 
ology  had  some  ability  to  simulate  the  effects  of  ice  interaction.  How- 
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ever,  a  particular  defect  of  this  rheology  Is  that  with  only  a  shear  vis¬ 
cosity  It  has  no  resistance  to  convergence  (see,  e.g. ,  Rothrock  1975a). 


Empirical  studies  using  a  more  "general*'  linear  viscous  law  containing 
both  bulk  and  shear  viscosities  (and  hence  having  compressive  strength,  see 
eq  9-11)  have  successfully  simulated  observed  ice  drift  and  deformation  far 
from  shore,  both  over  short  time  intervals  (Hibler  1974)  and  over  a  season¬ 
al  cycle  (Hibler  and  Tucker  1979).  In  particular,  Hibler  (1974)  showed 
that  a  linear  viscous  rheology  could  explain  the  observation  that  the  ice 
pack  often  tends  to  converge  in  a  low  pressure  system  in  the  winter  and 
diverge  in  the  summer  under  similar  forcing.  Best  es^mates  ^f  the  viscos¬ 
ity  magnitudes  yielded  viscosities  of  the  order  of  10 1  kg  s“  ,  with  the 
bulk  viscosity  (c)  about  twice  as  large  as  the  shear  viscosity  (n).  In  a 
seasonal  study  Hibler  and  Tucker  (1979)  showed  the  "best  fit"  viscous  pa¬ 
rameters  to  vary  in  a  regular  seasonal  manner.  This  seasonal  variation  was 
commensurate  with  the  physical  notion  that  sea  ice  has  less  thin  ice  in 
winter  (due  to  the  higher  freezing  rates)  and  hence  has  higher  strength. 

The  seasonal  viscosity  also  tended  to  correctly  simulate  observed  seasonal 
variations  of  the  angle  of  ice  drift  relative  to  the  wind. 

However,  while  generally  performing  well  far  from  shore,  the  linear 
viscous  rheology  has  one  major  problem  in  that  viscosities  for  reasonable 
simulations  in  the  central  pack  are  quite  different  from  viscosities  needed 
to  simulate  near-shore  behavior  (Rothrock  1975a).  For  example,  viscosities 
(two-dimensional)  of  the  order  of  1011  to  1012  kg  s-1  (Hibler  1974,  Hibler 
and  Tucker  1979)  are  needed  to  model  pack  ice  behavior  far  from  shore. 

Ne^ir  shore,  however,  the  effective  viscosities  may  be  as  small  as  108  kg 
s_1  (Hibler  ^t  al.  1974a),  and  in  very  small  channels  they  can  be  as  small 
as  105  kg  s  (Sodhi  and  Hibler  1980).  Such  variations  suggest  that  inclu¬ 
sion  of  some  type  of  nonlinear  behavior  is  necessary. 

An  additional  problem  with  linear  viscous  laws  is  that  their  physical 
basis  is  questionable.  In  particular,  the  local  characteristics  of  sea  ice 
appear  to  be  nonlinear  and  plastic  in  nature.  However,  even  though  sea  ice 
may  not  exhibit  linear  viscous  characterisics  on  a  small  scale,  Nye  (1973) 
has  noted  that  such  linear  laws  still  might  be  appropriate  for  a  time-  or 
spatially  averaged  response.  A  particular  case  of  this  concept  has  been 
proposed  by  Hibler  (1977),  who  suggested  that  on  appropriate  time  scales 
stochastic  fluctuations  in  the  deformation  rate  might  cause  the  average 
stress-strain  rate  relation  to  be  essentially  linear.  This  concept  may 
have  application  to  the  marginal  ice  zone  and  is  discussed  in  more  detail 
later  in  this  chapter. 

To  avoid  the  inherent  deficiencies  in  linear  viscous  rheologies  and  to 
better  describe  the  local  real-time  behavior  of  sea  ice,  Coon  (1974)  and 
Rothrock  (1975a)  proposed  a  plastic  rheology  for  sea  ice.  As  discussed 
earlier,  in  a  plastic  constitutive  law  the  stress  state  of  the  ice  is  pre¬ 
sumed  to  be  independent  of  the  magnitudes  of  the  stress  and  strain  rate. 
Such  a  law  allows  highly  nonlinear  behavior,  which  is  helpful  in  explaining 
ice  flow  both  near  and  far  from  shore.  In  addition,  it  contains  a  simple 
way  to  specify  a  low  tensile  strength  concurrently  wi th  a  high  compressive 
strength.  To  justify  such  a  plastic  law  physically,  Coon  and  Rothrock  ar¬ 
gued  that  such  behavior  is  consistent  with  the  local  physics  of  sea  ice. 

The  essential  argument  for  rate  Independence  is  that  the  work  done  in  ice 
deformation  is  primarily  due  to  ridging.  The  energy  expended  in  ridging, 


10 


however,  would  be  expected  to  be  independent  of  the  rate  of  ridge  building. 
This  concept  can  be  extended  into  a  complete  explanation  (in  terms  of  ridg¬ 
ing  and  frictional  losses)  of  all  energy  dissipated  during  plastic  flow. 
Thi6  particular  approach  is  discussed  in  the  next  section.  With  regard  to 
small  tensile  strengths,  it  was  argued  that  the  presence  (over  a  large 
scale)  of  many  cracks  prevented  any  effective  cohesion  under  dilation. 

While  the  concept  of  an  ideal  rigid  plastic  law  is  useful  theoretical¬ 
ly,  it  is  difficult  to  implement  in  practice.  In  particular,  to  examine  a 
plastic  rheology  numerically,  some  approximation  mist  be  made  to  describe 
the  “rigid"  behavior,  i.e.  the  behavior  when  the  ice  is  not  flowing  plasti¬ 
cally.  Two  approaches  have  been  suggested:  the  elastic-plastic  approach 
(Pritchard  1975)  and  the  viscous-plastic  approach  (Hibler  1979).  In  the 
elastic-plastic  method  the  ice  is  assumed  to  behave  elastically  for  certain 
strain  states.  The  advantage  of  this  is  that  it  allows  rigid  behavior  in 
the  sense  that  the  ice  can  sustain  a  variety  of  stresses  while  possibly 
undergoing  no  deformation.  However,  inclusion  of  the  elasticity  is  inher¬ 
ently  unwieldy  and  has  the  disadvantage  of  adding  elastic  waves  to  the  be¬ 
havior  of  the  model  (Colony  and  Rothrock  1980). 

In  the  viscous-plastic  approach  (Hibler  1979)  ,  sea  ice  is  considered 
to  be  a  nonlinear  viscous  continuum  characterized  by  both  nonlinear  bulk 
and  shear  viscosities  and  a  pressure  term.  These  nonlinear  viscosities  are 
adjusted  to  give  plastic  flow  for  normal  strain  rates.  To  close  the  sys¬ 
tem,  for  very  small  deformation  rates  the  ice  is  assumed  to  behave  as  a 
very  stiff  linear  viscous  fluid.  Thus  the  rigid  behavior  is  approximated 
by  a  state  of  very  slow  flow.  This  technique  is  computationally  advanta¬ 
geous.  Its  disadvantage,  however,  is  that  ice  that  is  actually  stationary 
may  be  modeled  as  slowly  creeping. 

Simulations  using  plastic  rheologies  have  reproduced  many  aspects  of 
sea  ice  dynamics.  Near-shore  simulations  by  Pritchard  et  al.  (1977)  have 
shown  the  capability  of  the  elastic-plastic  rheology  to  model  both  rela¬ 
tively  stationary  behavior  under  stress  and  intense  shearing  behavior  near 
coasts.  With  respect  to  seasonal  simulations,  Hibler  (1979)  has  shown  that 
in  addition  to  producing  near-shore  shear  zone  effects,  a  plastic  rheology 
yields  geographical  ice  thickness  buildup  and  ice  outflow  from  the  Arctic 
Basin  which  agree  well  with  observations.  Also,  for  sufficiently  high 
strengths,  Hibler  (1980a)  found  the  modeled  ice  cover  in  the  Arctic  Basin 
stopped  moving,  even  though  the  wind  forcing  is  significant.  Such  a  state 
of  no  motion  in  a  constricted  region,  even  though  the  external  forcing  is 
significant,  is  an  important  manifestation  of  a  nonlinear  rheology.  This 
characteristic  is  a  salient  feature  of  ice  drift  in  regions  with  narrow 
land  constrictions  such  as  the  Bay  of  Bothnia  or  the  Great  Lakes  of  the 
United  States. 

Analytic  analysis  of  a  one-dimensional  system 

To  get  some  feeling  for  the  relative  behavior  of  linear  versus  non¬ 
linear  rheologies  and  to  examine  the  effects  of  nonlinear  ice  interaction 
on  ice  drift  rates  it  is  useful  to  analyze  a  special  case  of  the  momentum 
balance.  Consider  the  one-dimensional  case  of  the  momentum  balance  employ¬ 
ing  only  a  linear  water  drag  term  cu,  an  external  constant  wind  stress  t, 
and  a  one-dimensional  ice  interaction  stress  o: 
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Consider  the  region  bounded  by  rigid  walls  at  x  ■  0  and  x  ■  L.  Let  us  now 
examine  the  linear  viscous  case  where  a  -  ?  3u/3x  with  €  a  constant.  Sub¬ 
stituting  this  expression  for  o  Into  eq  17  gives 


cu  -  C 


A  general  solution  of  eq  18  Is 

u  -  A  exp[(c/ c)^2x]  +  B  exp  [- (c/ c)1^2x  ]  +  ^  .  (19) 

Utilizing  the  boundary  conditions  u(0)  -  u(L)  *  0  we  obtain 

U(X)  *  2cslnh( AL)  Ue_AL  “  1  1  +  t1  -  eALl  e_AXi  +  i  (20) 

where  A  ■  (c/c)*^. 

Let  us  now  consider  a  rigid  plastic  case  with  constant  strength.  For 
our  plastic  rheology  we  will  assume  that  a  -  -P  for  3u/ 3x  <  0  and  o  ■  0  for 
3u/3x  >  0.  These  assumptions  define  a  rigid  plastic  rheology  with  no 
tensile  strength.  A  solution  for  this  case  may  be  constructed  by  noticing 
a)  that  for  any  deformation  o  -  -P  while  for  no  motion  0  a  >_  - P,  and  b) 
that  the  maximum  force  would  be  expected  to  take  place  at  the  right-hand 
boundary  since  the  wind  stress  has  built  up  at  this  point.  With  these 
considerations  in  mind  we  see  first  that  if  t  L  <  P  no  motion  of  any  kind 
will  occur  and  the  system  will  be  rigid.  If,  on  the  other  hand,  x  L  >  P 
then  a  solution  of  eq  17  can  occur  for  a  o  satisfying  the  above  plastic 
assumptions  and  o(L)  ■  -P.  Integrating  eq  17  from  x  to  L  we  have 

(cu  -  x)  -  -P  -  o(x)  .  (21) 

However,  we  must  also  have  3u/ 3x  *  0  except  at  the  boundary,  otherwise  o(x) 
*  -P  and  eq  17  cannot  be  satisfied.  Therefore,  u  is  a  constant,  as  is  x, 
and  we  obtain 

o(x)  -  -P  +  ( cuq  -  x)  (x  -  L).  (22) 

The  value  of  Uq  may  be  obtained  by  noting  that  the  total  force  acting  on 
the  rigid  block  moving  is  -cu0L  +  xL  -  P;  since  there  is  no  acceleration 

%  -  (23) 


Another  way  to  see  this  is  to  note  that  a  (x*0)  ■  0  since  divergence  is 
occurring  there  and  we  have  assimed  no  tensile  strength.  With  a  little 
thought,  this  type  of  simple  analysis  can  be  extended  to  a)  the  equations 
with  inertial  terms  (rigid  block  motion  occurs  with  the  block  gradually 
speeding  up  to  equllibrlun  with  a  spin-up  time  based  on  the  water  drag). 


and  b)  a  case  with  strong  wind  gradients  (for  sufficiently  large  gradients 
the  region  of  convergence  will  occur  over  a  finite  distance  rather  than  at 
a  single  location) .  In  this  latter  case  once  convergence  has  occurred  the 
stress  will  be  constant,  and  hence  free  drift  will  occur  even  though  the 
ice  Is  interacting.  These  more  general  cases  will  be  illustrated  with  some 
numerical  examples  discussed  below. 

This  analysis  is  also  easily  extended  to  include  nonlinear  water  drag 
terms,  an  extension  that  makes  it  more  comparable  to  the  free  drift  analy¬ 
sis  carried  out  earlier  in  the  momentin  balance  section.  Using  the  nota¬ 
tion  of  eq  6,  the  steady-state  momentum  balance  for  this  one-dimensional 
system  is 

paca^g  “  Pw°w^  “  ■3^*  (24) 

By  the  same  methods  as  employed  in  the  linear  case  for  a  constant  wind  the 
solution  for  this  system  is 

/ pacau!  “  P/L 

U  -/  a  8  -  (25) 

/  pwcw 

2 

if  P  <  Pa  ca  Ug  L,  and  U  *  0  otherwise.  For3a  numerical  comparison  rele¬ 
vant  to  the  Arctic  Basin  we  take  L  »  2.5  *  l(j  km,  which  is  the  scale  of 
the  Arctic  Basin,  and  P  ■  8.25  *  lO4  N  a"1,  an  empirical  best  fit  ice 
strength  for  3-m-thick  ice  in  the  Arctic  Basin  in  winter  (Hibler  and  Walsh 
1982).  Using  these  numbers  together  with  the  drag  coefficients  used  in 
Figure  3,  we  obtain  a  comparison  of  free  drift  and  drift  with  Ice  interac¬ 
tion  as  shown  in  Figure  6.  Note  that  for  large  wind  speeds  there  is  little 


Figure  6.  Comparison  of  ice  drift  rates 
for  an  idealized  one-dimensional  system 
with  and  without  plastic  ice  interaction. 
Quadratic  boundary  layers  were  used,  and 
for  the  ice  interaction  case  typical 
strength  and  length  scales  relevant  to 
the  Arctic  Basin  were  assumed.  Other 
constants  are  the  same  as  used  in  Figure 
3  for  3-m-thick  ice. 


difference  between  free  drift  and  the  ice  interaction  drift  whereas  for 
smaller  wind  speeds  the  difference  is  very  marked,  with  the  rigid  plastic 
system  totally  stopping.  A  significant  difference  from  the  linear  water 
drag  case  is  that  the  slope  of  the  ice  drift  rate  versus  wind  speed  de¬ 
creases  asymptotically  to  the  free  drift  slope  for  large  geostrophic  wind 
speeds.  In  the  one-dimensional  linear  water  drag  case  the  free  drift  and 
ice  interaction  curves  would  be  parallel  but  would  have  different  x-axis 
intercepts. 

A  mirror  reflection  of  this  figure  holds  for  negative  ice  velocities. 
Consequently,  if  one  were  to  force  this  model  with  a  zero  mean  Gaussian 
distribution  of  wind  speeds  and  then  fit  the  results  with  a  linear  model, 
the  effect  of  the  ice  interaction  would  be  to  reduce  the  ratio  of  ice  drift 
to  wind  speed.  This  ratio  would,  however,  be  dependent  on  the  variance  of 
the  wind  speeds,  with  smaller  variances  yielding  more  of  a  slope  reduction. 
These  considerations  are  relevant  to  linear  correlations  of  ice  drift  rates 
with  winds  (Thondike  and  Colony  1982),  which  yield  only  slight  reductions 
in  winter.  This  simple  analysis  suggests  that  part  of  the  reason  for  such 
results  may  be  the  increased  intensity  and  hence  variance  of  winds  in  win¬ 
ter.  This  analysis  also  Indicates  that  although  useful,  such  linear  corre¬ 
lations  may  be  masking  some  of  the  essential  physics  of  ice  drift. 

Some  numerical  solutions  for  an  idealized  geometry 

To  examine  the  behavior  of  a  plastic  rheology  for  a  more  realistic 
two-dimensional  geometry  and  for  spatially  varying  wind  forcing  some  numer¬ 
ical  solutions  were  carried  out  using  the  viscous-plastic  approach  to  mod¬ 
eling  plastic  flow.  These  solutions  help  verify  some  of  the  simple  con¬ 
cepts  mentioned  above  in  the  one-dimensional  analysis  and  also  demonstrate 
the  utility  of  the  viscous-plastic  approach  for  modeling  plastic  flow. 

For  the  momentum  balance  eq  1  is  used  with  the  geostrophic  current 
Uw  set  equal  to  zero  and  with  the  momentum  advection  term  set  equal  to 
zero: 

3u 

m  x  u  +  Ta  +  tm  +  F.  (26) 

For  modeling  the  ice  interaction  the  ice  is  considered  to  have  a  nonlinear 
viscous  plastic  constitutive  law  similar  to  that  mentioned  earlier. 

°ij  “  2T1<Vj.P)  z±i  +  UC^.P)  -  n  C  ,P) ]  efck  6^  -  P o1  j / 2  (27) 


where  oy  is  the  two-dimensional  stress  tensor,  the  strain  rate  tensor, 
tensor,  P*/2  a  pressure  term,  and  ?  and  n  are  nonlinear  bulk  and  shear  vis¬ 
cosities.  For  normal  strain  rates  the  dependences  of  C  and  n  on  and  P* 
is  given  by  eq  13-16  so  that  the  stress  state  lies  on  an  elliptical  yield 
curve  passing  through  the  origin  with  a  no-stress  condition  applying  for 
pure  divergence.  For  very  small  strain  rates  the  viscosities  specified  by 
eq  13-16  become  arbitrarily  large.  To  avoid  this  they  are  chosen  to  be  the 
minimum  of  the  plastic  values  and  some  large  limiting  values  dependent  on 
the  ice  strength  P*.  For  the  calculations  performed  here  the  standard  lim¬ 
iting  values  were  taken  to  be 


max 


C  led 
max 


(29) 


Sensitivity  tests  (Hibler  et.  al.  1981)  show  these  simulation  results  to  be 
relatively  insensitive  to  these  maximian  viscosities  unless  they  are  taken 
to  be  about  two  orders  of  magnitude  smaller.  However,  the  maximum  viscosi¬ 
ty  values  do  affect  computational  speeds,  especially  at  finer  resolution 
(see  Hibler  et  al.  1981). 

To  solve  eq  26-29  and  12-16  numerically,  a  time  stepping  procedure 
using  finite  difference  techniques  is  employed.  The  computer  code  for  this 
procedure  was  documented  by  Hibler  (1980b).  Briefly,  the  numerical  proce¬ 
dure  works  as  follows.  The  viscosities  used  in  the  momentum  balance  are 
based  on  the  deformation  field  from  the  previous  time  step.  Using  these 
viscosities  a  new  velocity  field  is  obtained,  a  new  set  of  viscosities  are 
estimated,  and  another  linearized  equation  solved.  We  will  refer  to  each 
of  these  relaxation  solutions  as  an  "iterative  time  step."  By  carrying  out 
several  iterative  time  steps  at  each  "physical  time  step,"  ideal  plastic 
flow  may  be  approached.  In  the  simulations  in  this  section  15  iterative 
time  steps  are  carried  out  at  each  physical  time  step,  so  that  plastic 
equilibrium  is  approached  at  each  time  step.  The  effect  of  smaller  numbers 
of  iterative  time  steps  is  discussed  in  more  detail  by  Hibler  et  al. 

(1981). 


Using  these  equations  and  numerical  method  a  series  of  simulations 
over  a  region  148  by  148  km  with  fixed  boundaries  and  fixed  geos  trophic 
wind  in  the  x-direction  were  carried  out  (see  Fig.  7).  The  boundary  condi¬ 
tions  consisted  of  specifying  zero  velocity  at  the  boundaries.  For  later 
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Figure  8.  Equilibrium  velocity  field  for 
the  standard  simulation. 


Table  1.  Numerical  parameters  used  in  standard  case  simulation. 

Ca  *  0.0012  pa  *  1.3  kg  n'3  At  »  3  hours  Ug  *»  8  m  s"1 

C w  =»  0.005  P*  =  5000  N  m_1  Cnax  =  (5*10?s)P*  Vg  =•  0.0 

e  =  2  m  =>  0.91- 10 3  kg  m-2  n^x  »  ^ma x/e2  uw  “  0.0 

f  -  1.46*10~4s-1  Ax  -  Ay  =  18.5  km  4>  =  0  »  25°  Vw  -  0.0 


reference  the  standard  simulation  utilized  the  parameters  shown  in  Table  1. 
The  equilibrium  velocity  field  for  the  standard  simulation  is  shown  in 
Figure  8.  While  there  are  some  deviations  due  to  the  two-dimensional  char¬ 
acter  of  the  rheology  and  boundaries  this  result  basically  represents  a 
rigid  block  motion  that  is  consistent  with  the  simple  one-dimensional  argu¬ 
ments  presented  in  the  previous  subsection.  In  this  block  motion  the  main 
converging  deformation  occurs  at  the  right-hand  boundary,  while  divergence 
occurs  at  the  left-hand  boundary.  There  is  also  some  slight  reduction  of 
velocity  at  the  side  boundaries  due  to  shear  stresses. 

Note  that,  as  discussed  above,  the  plastic  solution  is  qualitatively 
very  similar  to  free  drift  (shown  in  Figure  8  as  a  single  vector  ouside  the 
grid)  in  that  the  whole  sea  ice  region  moves  practically  at  the  same  veloc¬ 
ity.  However,  there  are  quantitative  differences.  As  in  the  one-dimen¬ 
sional  case  (see  Fig.  6),  a  major  difference  is  a  smaller  velocity  magni¬ 
tude  in  the  plastic  system  than  for  free  drift.  A  second  significant  dif- 
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a.  The  component  of  the  enter  point  velocity  versus 
time. 


Distance  (km) 


b.  Profiles  (taken  parallel  to  the  x  axis  and 
through  the  center  point)  of  the  x  velocity 
componen  t. 

Figure  9.  Inertial  spin-up  characteristics  of  idealized 
viscous-plastic  sea  ice  model. 


ference  is  a  substantially  reduced  turning  angle  (  15°  reduction)  between 
the  plastic  ice  velocity  vectors  and  the  surface  winds.  As  discussed  in 
the  earlier  sections  on  the  momentum  balance  (see  Fig.  4)  this  can  be  ex¬ 
plained  by  the  fact  that  the  force  due  to  internal  ice  stress  approximately 
opposes  the  ice  motion.  Moreover,  as  the  ice  strength  increases  this  turn¬ 
ing  angle  difference  becomes  even  more  pronounced. 

The  spin-up  results  for  the  standard  case  are  shown  in  Figure  9.  In 
these  spin-up  results  3-minute  physical  time  steps  were  used  as  compared  to 
3  hours  used  in  the  other  tests.  Figure  9a  shows  the  x  component  of  the 
center  point  velocity  versus  time,  while  Figure  9b  shews  velocity  profiles 
(along  the  x  axis)  at  9-minute  intervals.  As  can  be  seen  the  spin-up  is 
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Figure  10.  Equilibrium  x-veloci ty 
component  of  the  center  point  as  a 
function  of  ice  strength  and  for 
different  yield  curves.  Curves 
are  shown  for  the  standard  yield 
curve  (e=2)  and  for  a  narrower 
elliptical  yield  curve  with  less 
shear  strength  (e=4). 


also  a  rigid  motion,  with  the  block  motion  gradually  increasing  in  accor¬ 
dance  with  the  magnitude  of  the  inertial  term.  A  notable  feature  of  these 
numerical  spin-up  results  is  that  they  are  free  of  elastic  wave  effects 
that  result  from  elastic-plastic  rheologies  (see,  e.g.,  Colony  and  Rothrock 
1980).  As  a  result,  the  viscous  plastic  approach  allows  a  more  direct  anal¬ 
ysis  of  the  essential  physics  of  ice  interaction  with  minimal  numerical 
artifacts.  Tests  were  also  done  with  different  strengths.  Spin-up  times 
for  these  cases  were  slightly  different,  due  to  the  fact  that  different 
equilibrium  velocities  resulted  in  different  water  stress  terms.  However, 
with  a  linear  water  drag  the  spin-up  times  were  essentially  the  same  for 
different  strengths. 

The  response  of  this  simple  system  was  also  examined  for  a  variety  of 
forces  of  different  magnitudes  and  having  spatial  variation.  Analysis  of 
the  results  shows  that  the  equilibriun  velocity  characteristics  depend  on 
the  physical  size  of  the  basin,  external  forcing,  and  ice  strength.  From 
dimensional  analysis  it  can  be  deduced  that  near  the  "motion"  to  "no-mo- 
tion"  transition  a  key  dimensionless  parameter  is  P*/(Ta*L)  =  B,  where  L 
is  the  size  of  the  region,  P*  the  ice  strength  and  Ta  the  wind  stress. 

To  examine  the  dependence  on  this  parameter  a  series  of  tests  were  done 
with  fixed  resolution  (Ax),  size  of  area  (L),  and  wind  stress  ( Ta) .  In 
each  test  a  different  value  of  P*  (fixed  in  time  and  space)  was  used.  All 
other  parameters  were  the  same  as  the  standard  case.  To  verify  the  scal¬ 
ing,  tests  were  also  done  with  L  and  xa  varying,  with  similar  results  as 
long  as  8  was  fixed.  (To  obtain  identical  results  the  water  drag  and  ice 
mass  would  also  have  to  be  modified.)  Figure  10  shows  the  component  of  the 
center  point  velocity  versus  strength  from  these  simulations.  This  figure 
shows  the  strength  for  which  the  system  behaves  rigidly  and  also  shows  that 
once  the  strength  is  low  enough  to  allow  the  ice  to  move,  halving  the 
strength  brings  the  velocity  field  almost  to  free  drift  values.  The  value 
of  B  for  the  rigid  "motion"  to  "no-motion"  transition  is  about  half  of  the 
value  estimated  for  the  one-dimensional  case.  This  reduced  value  of  B  is 
primarily  due  to  the  effect  of  shear  stress,  which  was  not  considered  in 
the  one-dimensional  analysis.  This  hypothesis  is  verified  in  Figure  10, 
which  shows  equilibrium  velocities  versus  P*  for  a  smaller  ratio  of  shear 
to  compressive  strength  (e  =  4). 
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Figure  11  shows  the  effect  of 
a  spatially  varying  wind  field. 
Illustrated  in  part  a  of  this  fig¬ 
ure  is  the  equilibrium  ice  velocity 
field  for  a  wind  field  varying  di¬ 
rectly  from  16  m  s“ 1  at  the  left- 
hand  boundary  to  zero  at  the  right- 
hand  boundary.  As  can  be  seen, 
this  large  velocity  field  gradient 
causes  the  region  of  rigid  motion 
to  recede  to  the  left-hand  side, 
with  the  convergence  now  occurring 
over  a  larger  region  at  the  right- 
hand  side.  This  convergence  has  to 
occur  over  a  band  because  the  gra¬ 
dient  of  the  stress  in  the  rigid 
region  now  must  change  in  order  to 
match  the  external  force,  and  the 
changes  available  for  a  given  ice 
strength  are  limited.  Figure  lib 
shows  the  case  where  there  is  an 
offshore  wind  gradient  with  wind 
varying  from  zero  the  left-hand 
boundary  to  16  m  s~  at  the  right- 
hand  boundary.  In  this  case  the 
region  of  divergence  is  enlarged. 


Figure  11.  Effect  of  gradients  in 
the  wind  fields  on  the  equilibrium 
velocity  field  for  the  idealized 
domain  in  Figure  7.  In  part  (a)  the 
wind  is  in  the  positive  x  direction 
and  varies  from  16  m  s"  at  the  left- 
hand  boundary  to  0  at  the  right-hand 
boundary.  In  part  (b)  the  wind  (also 
in  the  positive  x  direction)  starts 
at  0  at  the  left-hand  boundary  and 
goes  to  16  m  s“  at  the  right-hand 
boundary. 


An  interesting  feature  of  Fig¬ 
ure  11  is  the  symmetry  between  the 
two  ice  drift  curves  for  opposite 
gradients  in  the  wind  fields.  This 
occurs  even  though  the  diverging 

ice  has  effectively  no  stress,  whereas  the  converging  ice  has  a  stress 
state  at  the  plastic  yield.  However,  in  the  converging  case  even  though 
the  stress  is  large  there  is  effectively  no  stress  gradient  and  hence  no 
significant  force  due  to  ice  interaction.  In  general  this  figure  demon¬ 
strates  some  of  the  counter-intuitive  results  that  can  occur  with  a  non¬ 
linear  ice  interaction.  In  this  case,  even  though  the  ice  is  interacting 
strongly  the  local  drift  still  obeys  free  drift. 


These  gradient  results,  in  conjunction  with  the  spin-up  results,  have 
application  to  rapidly  moving  fronts.  Such  fronts  contain  large  wind  gra¬ 
dients,  and  on  the  basis  of  these  simulations  would  be  expected  to  cause 
moving  zones  of  deformation  to  propagate  ahead  of  the  front.  The  block 
motion  results  also  give  insights  into  why  simple  free  ice  drift  ice  fore¬ 
casts  have  generally  been  quite  useful  (see,  for  example,  Zubov  1943). 
Basically,  in  the  rigid  motion  case  the  rigid  block  behaves  very  much  like 
freely  drifting  ice,  except  that  the  plastic  stress  tends  to  renormalize 
the  wind  field  (or  whatever  external  forcing  exists)  to  a  smaller  value. 
Also,  such  free  drift  rates  are  usually  invoked  to  predict  the  drift  and 
not  the  deformation,  and  hence  are  primarily  based  on  the  smoothly  varying 
components  of  the  wind  fields. 


Figure  12.  Simulated  and  observed  trajec¬ 
tories  of  a  buoy  off  East  Greenland,  assum¬ 
ing  either  a  viscous  plastic  ice  interac¬ 
tion  or  no  ice  interaction. 


Stochastic  effects  and  an  unconstrained  pressure  term 

Far  from  the  ice  edge,  where  the  floes  are  large  and  closely  packed, 
the  ideal  rigid  approximation  is  probably  good.  However,  near  the  edge 
there  may  be  a  great  deal  of  random  bumping  of  floes.  This  in  turn  may  in¬ 
troduce  an  unconstrained  pressure  terra  which  can  do  work  on  the  continuum 
system.  While  such  a  rheology  is  still  a  special  case  of  the  viscous  plas¬ 
tic  rheology,  it  can  cause  a  different  type  of  ice  edge  dynamics  than  is 
yielded  by  rigid  plastic  analysis.  However,  random  bumping  is  likely  lo¬ 
calized  at  the  edge  and  will  be  superimposed  on  a  larger  scale  behavior 
that  has  more  typical  plastic  characteristics.  Specifically,  on  the  large 
scale  it  is  likely  that  the  ice  interaction  has  a  rectifying  effect  on  the 
ice  motion.  Under  off-ice  winds,  for  example,  the  plastic  nature  of  the 
ice  will  cause  it  to  diverge  relatively  freely,  whereas  under  on-ice  winds 
the  ice  interaction  will  prevent  further  convergence. 

Preliminary  simulations  of  the  East  Greenland  pack  (Tucker  and  Hibler 
1981)  have  demonstrated  such  a  rectifying  effect,  and  indicate  that  it  is 
an  important  factor  in  the  tendency  of  the  ice  to  drift  parallel  to  the 
coast.  Figure  12,  for  example,  shows  simulated  and  observed  trajectories 
for  a  buoy  off  East  Greenland,  assuming  a  viscous-plas tic  ice  interaction. 
As  can  be  seen,  the  ice  interaction  tends  to  prevent  onshore  motion  of  the 
buoy,  even  though  the  winds  tend  to  drive  it  against  the  coast. 

One  way  to  examine  the  effect  of  random  bunping  and  other  stochastic 
fluctuations  is  to  assume  a  plastic  law  locally  and  then  average  over  sto¬ 
chastic  effects  in  space  and  time.  Such  an  approach  was  employed  by  Hibler 
(1977)  in  an  attempt  to  relate  plastic  law  parameters  to  those  used  in  vis¬ 
cous  models.  To  demonstrate  the  general  idea  let  us  consider  the  specific 
case  of  a  rigid  plastic  law  having  a  circular  yield  curve.  Referring  back 
to  eq  9-16  the  stress  strain  rate  relation  for  such  a  system  is: 


Oy  -  <P*/A)  [e  ]  -  6  P*/2 


where 


A  =  [(e2  +  e2  )  +  4e2  +  2e  t  )1/2.  (31) 
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To  obtain  a  statistical  average,  assume  that  ey  fluctuates  stochas¬ 
tically  about  a  mean  value  (the  fluctuations  could  be  in  either  time 

or  space)  with  some  Gaussian  probability  density 

G( Ey )  =  [X(2ir)l/2]-1  exp  [-<e  -  <e  »2/2X2l  (32) 

where  X  is  the  standard  deviation  of  the  strain  rate  fluctuations.  Conse¬ 
quently,  the  statistical  average  of  Oy ,  denoted  by  <Ojj>,  is  given  by 

<°ij>  =  t  A3(27r)3/2]-1  /  exp  {-l(en  -  <en»2  +  (  e22  -  <e22»2 

+  ( £ i 2  "  <C12>)^ 1/2X2 }o  (ei i ,£22,el2)d£i ide22d£i2  (33) 

where  Ojj  (£11,^22.^12)  is  given  by  eq  30.  From  this  equation,  two  limit¬ 
ing  cases  are  clear.  In  particular,  for  very  small  X,  which  corresponds  to 
very  small  fluctuations  in  f-jj,  EjH  in  may  be  replaced  by  <5^j>,  in 
which  case  the  plastic  law  is  obtained  for  the  average  values.  For  very 
large  X,  which  corresponds  to  large  fluctuations  compared  to  the  mean  val¬ 
ue,  the  right-hand  side  of  eq  33  may  be  expanded  in  powers  of  <ejj>/X,  a 
procedure  which  would  yield  to  first  order  a  viscous-like  law,  i.e. 
proportional  to  <£ij>.  It  is  worth  noting  here  that  this  expansion  argu¬ 
ment  is  independent  of  the  type  of  distribution  and  constitutive  law  as¬ 
sumed.  For  example,  considering  an  arbitrary  constitutive  law  and  a  gener¬ 
al  probability  distribution,  we  could  still  expand  in  powers  of  <£jj>/X, 
although  the  actual  coefficients  would  depend  on  the  functional  law  as¬ 
sumed. 


For  the  special  case  of  a  circular  yield  curve  it  is  possible,  by 
using  eq  33,  to  estimate  analytically  the  coefficients  of  the  linear  terms 
in  <£jLj>/X.  Substituting  eq  30  and  31  into  eq  33  and  changing  variables, 
after  some  algebra  (see  Hibler  1977),  one  obtains 


<n^>  -  P*  {0 .35  Ce^VX  -  0.5o  }  (34) 

where  0.35  is  an  approximation  to  two  significant  figures  of  a  combination 
of  logarithms  of  integers  and  pi.  Referring  back  to  eq  9 ,  this  expression 
is  basically  a  linear  viscous  law  with  a  pressure  term  and  with  n  -  1 ;  *» 
0.173  P*/X. 


For  the  elliptical  yield  curve,  analytic  estimation  of  n  and  C  is  dif¬ 
ficult,  although  their  ratio  can  be  deduced.  In  particular,  it  is  clear  by 
symmetry  that  in  expanding  eq  33  in  powers  of  <e^j>/X  the  lowest-order 
terms  will  yield  an  average  stress-strain  relationship  of 
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(35) 


i 


<<7lj>  =  2TV(alj>  +  (C  -  Tl)<ekk>  6lj  “  6ij/2p*  (35) 

2 

where  rj n  »  e  .  Consequently,  as  the  ellipse  becomes  more  elongated,  the 
ratio  5/n  will  increase. 

For  a  more  specific  application  of  this  concept  one  can  consider  the 
case  of  fluctuation  in  time  of  the  deformation  at  a  given  location.  Using 
observed  and  simulated  deformation  results,  Hibler  (1977)  was  able  to  show 
that  averaging  over  a  period  of  only  a  few  days  yielded  a  more  linear-like 
behavior,  together  with  a  finite  pressure  term.  Considering  other  effects, 
such  as  inertial  oscillations  that  are  probably  prominent  near  the  edge,  it 
is  likely  that  necessary  averaging  intervals  for  these  stochastic  arguments 
will  be  quite  small,  and  likely  shorter  than  typical  temporal  wind  varia¬ 
tions. 

The  ramifications  of  such  a  pressure  term  and  nonlinear  behavior  near 
the  edge  are  not  clear.  However,  one  interesting  case  has  been  investi¬ 
gated  by  Roed  and  O’Brien  (1981).  They  considered  the  case  of  an  uncon¬ 
strained  pressure  term  only,  without  shear  strength  of  any  kind,  and  with¬ 
out  water  drag.  When  this  rheology  was  coupled  to  a  Coriolis  force  and  to 
a  simple  continuity  equation  for  the  ice  compactness,  geostrophic  adjust¬ 
ment  tended  to  cause  a  jet-like  effect  to  appear  at  the  ice  edge.  While 
this  calculation  is  a  bit  unrealistic  it  is  suggestive  of  some  of  the  rhe¬ 
ology-induced  phenomena  that  may  occur  near  the  ice  edge. 

Relative  magnitudes  of  shear  and  compressive  strengths 


The  above  stochastic  arguments  give  some  justification  for  using  lin¬ 
ear  viscous  laws  in  certain  circumstances.  One  example  of  such  a  calcula¬ 
tion  is  an  idealized  infinite  boundary  calculation  performed  by  Hibler 
(1974).  The  main  reason  for  this  calculation  was  to  compare  observed  de¬ 
formation  models.  However,  the  calculations  also  give  some  insight  into 
the  relative  magnitudes  of  shear  and  compressive  strengths.  Specifically, 
these  idealized  calculations  yield  a  best  fit  bulk  viscosity  about  twice  as 
large  as  the  shear  viscosity.  This  in  turn  suggests  that  although  of  the 
same  order  of  magnitude,  compressive  strengths  are  likely  larger  than  shear 
s  trengths. 

The  basic  idea  in  the  calculation  was  to  consider  a  linearized  momen¬ 
tum  balance  (linear  air  and  water  drag)  employing  a  linear  constitutive  law 


2TleiJ  +  <«  -  ">  PVk  V 


Solving  these  equations  in  the  case  of  an  infinite  boundary,  and  taking  the 
limit  of  large  viscosities,  yielded  the  result  that  the  divergence  rate 
over  the  vorticity  is  given  by 


TTTT)  (37) 

where  $  is  the  turning  angle  of  the  atmospheric  boundary  layer.  Comparison 
of  observed  ice  vorticity  and  divergence  rate  over  a  month-long  period  (see 
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Hibler  1974)  did  show  a  negative  correlation,  with  the  vortlclty  being 
about  twice  as  large  as  the  divergence  rate.  This  In  turn  Indicated  that  £ 
*  2n  and  suggests  that  the  compressive  strength  Is  larger  than  the  shear 
strength. 


RELATION  OF  ICE  STRENGTH  TO  ICE  THICKNESS  CHARACTERISTICS 

In  the  above  studies  we  have  discussed  the  behavior  of  Ice  dynamics 
under  conditions  where  a  constant  strength  Is  specified.  However,  in  prac¬ 
tice  the  ice  strength  is  related  to  ice  thickness  characteristics  in  some 
way.  To  approximate  the  concept  there  have  been  two  general  approaches  to 
date:  a  multi-level  approach  developed  by  Thorndike  et  al.  (1975)  and  a 
two-level  approach  developed  by  Hibler  (1979).  In  the  multi-level  approach 
the  ice  thickness  distribution  is  broken  down  into  a  large  number  of  levels 
and  the  strength  determined  by  a  hypothesized  model  of  the  ice  ridge  build¬ 
ing  process.  In  the  two-level  approach  the  ice  strength  is  empirically  re¬ 
lated  to  the  compactness  and  the  mean  ice  thickness.  In  the  next  two  sub¬ 
sections  we  discuss  these  approaches.  In  the  first  section  a  critique  of 
the  Thorndike  et  al.  (1975)  ridge  building  process  is  presented,  and  other 
redistributors  are  presented  and  discussed.  This  section  also  gives  esti¬ 
mates  of  ice  strength  based  on  mechanical  considerations  and  compares  these 
strengths  to  best  fit  strengths  used  in  numerical  simulations  to  obtain 
agreement  with  observed  drift.  The  two-level  formulation  of  ice  strength 
is  discussed  in  the  second  subsection.  A  brief  comparison  of  one-year  sim¬ 
ulations  of  the  Arctic  Basin  ice  cover  employing  the  multi-level  and  two- 
level  formulations  is  given  in  the  third  subsection.  Finally,  in  the 
fourth  subsection  some  of  the  small  scale  ramifications  of  a  nonlinear 
coupling  between  ice  strength  and  ice  thickness  characteristics  are  dis¬ 
cussed. 

Variable  thickness  parameterization  of  ice  strength 

Recent  work  on  the  strength  of  a  variable  thickness  ice  cover  has  been 
largely  based  on  the  concept,  first  suggested  by  Parmerter  and  Coon  (1972) 
in  a  kinematic  ridge  model,  that  the  major  work  done  in  ridge  building  is 
due  to  potential  energy  changes  caused  by  ice  pile-up.  Rothrock  (1975b) 
followed  up  this  idea  by  suggesting  that  this  concept  be  extended  to  two 
dimensions  by  insisting  that  two-dimensional  deformational  work  be  ex¬ 
plicitly  related  to  the  amount  of  work  done  by  ridging.  However,  if  one 
accepts  this  hypothesis  it  then  becomes  critical  to  determine  how  ice 
ridges. 

To  formulate  the  ridge  building  process  on  the  geophysical  scale 
Thorndike  et  al.  (1975)  suggested  that  thin  ice  be  redistributed  into 
thicker  ice  categories.  The  precise  manner  in  which  this  redistribution 
should  be  carried  out  is  not  clear.  As  an  initial  guess,  Thorndike  et  al. 
(1975)  suggested  a  redistributor  which  transfers  ice  into  categories  that 
are  fixed  multiples  of  the  initial  thickness.  Ridge  observations  and  theo¬ 
retical  considerations  suggest  that  such  a  redistributor  is  unrealistic. 

In  order  to  provide  a  more  realistic  parameterization  of  the  ridging  pro¬ 
cess,  Hibler  (1980a)  proposed  a  scaling  law  for  ridge  building,  and  used 
this  scaling  law  in  a  seasonal  equilibrium  simulation  of  the  Arctic  Basin 
ice  cover.  In  the  discussion  below  we  examine  in  more  detail  the  strength 
characteristics  of  different  redistribution  processes,  and  compare  these 
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processes  to  ridge  morphological  data.  In  addition,  annual  ridge  build-up 
results  from  the  seasonal  equilibrium  simulation  of  Hibler  (1980a)  are 
reported  and  compared  to  observations  of  ice  roughness. 

General  framework.  Following  Thorndike  et  al.  (1975)  the  ice  thick¬ 
ness  characteristics  are  described  by  an  areal  ice  thickness  distribution 
g(h,x,t),  where  g(h,x,t)dh  is  the  fraction  of  area  (in  a  region  centered  at 
position  x  at  time  t)  covered  by  ice  with  thickness  between  h  and  h  +  dh. 
This  distribution  evolves  in  response  to  deformation,  advection,  growth  and 
decay.  Neglecting  lateral  melting  effects  Thorndike  et  al.  (1975)  derived 
the  following  governing  equation  for  the  thickness  distribution: 

|f  +  V.(ug)  +  -  *  (38) 

where  f  is  the  vertical  growth  (or  decay)  rate  of  ice  of  thickness  h  and  ♦ 
is  a  redistribution  function  (depending  on  h  and  g)  that  describes  the  cre¬ 
ation  of  open  water  and  the  transfer  of  ice  from  one  thickness  to  another 
by  rafting  and  ridging.  In  general,  f  is  considered  to  be  a  function  of 
g(h),  and  <|>  a  function  of  both  g(h)  and  the  rate  of  deformation.  Except 
for  the  last  two  terms,  eq  38  is  a  normal  continuity  equation  for  g.  The 
last  term  on  the  left-hand  side  can  also  be  considered  a  continuity  re¬ 
quirement  in  thickness  space  since  it  represents  a  transfer  of  ice  from  one 
thickness  category  to  another  by  the  growth  rates.  An  important  feature  of 
the  Thorndike  et  al.  (1975)  theory  is  that  it  presents  an  "Eulerian”  de¬ 
scription  in  thickness  space.  In  particular,  growth  occurs  by  rearranging 
the  relative  areal  magnitudes  of  different  thickness  categories. 

This  equation  can  be  formally  generalized  (Hibler  1980a)  to  include 
lateral  melt  by  simply  adding  a  source  and  sink  term  FL(g,h),  such  that 

Q  FLdh  «  0.  (39) 

Equation  39  follows  from  the  fact  that,  by  definition,  lateral  melting  of 
ice  will  be  compensated  for  by  a  change  in  the  extent  of  open  water.  In 
addition,  Fl  >  0  for  h  ■  0  and  Fl  <  0  for  h  >  0.  These  conditions  reflect 
the  fact  that  lateral  melting  decreases  the  areal  extent  of  ice  while  in¬ 
creasing  open  water  extent.  Adding  this  lateral  growth  term  to  eq  38,  the 
more  general  thickness  distribution  equation  becomes 

If-  +  Mug)  +  -  Fl  +  <|>.  (40) 

Mechanical  redistributor.  Consider  now  a  two-dimensional  continuum 
with  ice  described  by  a  plastic  rheology  such  that  the  stress  state 
0i1^i1»p*)  is  a  function  of  the  strain  rate  and  a  strength  P*.  Then, 
following  Thorndike  et  al.  (1975)  and  Rothrock  (1975b),  the  mechanical 
redistribution  function  <|>  can  be  reasonably  postulated  to  satisfy  the  con¬ 
straints 


Jq  <|dh  -  V.  vi 


(41) 


■.*  \' 
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Q  hi|idh  -  0  (42) 

C  Q  h2<Kh)dh  -  e  (43) 

where  u  is  the  ice  velocity  field  and  C  a  constant.  Considering  only  the 
work  done  on  the  ice  by  gravitational  and  buoyancy  forces,  C  would  be  given 
by 

,  P_  -  P,  . 


with  the  density  of  ice,  the  density  of  water,  and  g  the  acceleration 
due  to  gravity.  Equation  41  follows  from  the  constraint  that  ij>  renormalizes 
the  g  distribution  to  unity  due  to  changes  in  area.  Equation  42  follows 
from  conservation  of  mass  and  basically  states  that  does  not  create  or 
destroy  ice  but  merely  changes  its  distribution.  (An  additional  assumption 
in  eq  42  is  that  the  ice  mass  is  related  in  a  fixed  manner  to  the  thick¬ 
ness.)  The  third  constraint  (Rothrock  1975b)  specifies  that  work  done  in 
building  ridges  is  equal  to  the  deformational  work.  An  important  feature 
of  this  constraint  is  that,  for  an  arbitrary  plastic  yield  curve,  it  forces 
some  ridging  to  occur  even  though  there  may  be  no  net  convergence.  In  its 
present  form  (with  C  ■  C^)  eq  43  ignores  frictional  losses  occurring  in 
ridging  and  considers  only  the  potential  energy  change  due  to  deformed  ice 
being  lifted  against  gravity  and  being  forced  down  against  buoyancy.  How¬ 
ever,  Rothrock  (1975b)  has  estimated  these  frictional  losses  to  be  the  same 
order  of  magnitude  as  the  potential  energy  changes.  Consequently,  as  a 
crude  approximation  for  frictional  losses,  eq  43  is  retained  here  with  a 
constant  C  equal  to  2Cb> 

A  redistribution  which  satisfies  these  constraints  is  (using  the  con¬ 
vention  that  repeated  subscripts  are  summed  over) 

*-  fiOOKP*)-1  Oyty  +  ekk]  +  (P*)”1  o^e^WrOi.g)  (45) 


»  ,  «  *1 
•  ■■.■■‘.--'.•I 


where  5(h)  is  the  dirac  delta  function  and 


wr(h,g) 


[-P(h)g(h)  +  J*  Y(h*  ,h)P(h’)g(h’)dh'] 
/~[p(h)g(h)  -  Jq  Y(h’  ,h)P(h’  )g(h’  )dh’  ]dh 


The  first  term  in  eq  45  specifies  the  amount  of  open  water  created,  while 
the  second  term  describes  the  transfer  of  thin  ice  to  thicker  categories  by 
ridging.  In  this  formalism  Y(hi,h2)dh2  can  be  thought  of  as  the  area  of 
ice  put  into  the  thickness  interval  [h2,h2  +  dh2]  when  a  unit  area  of  ice 
of  thickness  hi  is  used  up.  Basically  the  integral  over  Y  specifies  where 
the  ice  is  transferred  to  during  ridging.  The  function  P(h),  on  the  other 
hand,  is  some  probability  function  specifying  which  categories  are  being 
destroyed  by  ridging.  Written  in  this  form  i|>  automatically  satisfies  eq 
41.  To  satisfy  eq  42  and  43  leads  to  the  additional  constraints,  respec- 
t ive ly , 


isi 

*'»*•  iN'kV 

*  •  .*  "J*k  , 


S  /  f  s 


(47) 


fQ  Y(h' ,h)hdh  -  h' 

C  /£  Wrh2dh  =  P*.  (48) 

Equation  47  requires  that  transfer  of  ice  from  one  category  to  another  does 
not  change  its  mass,  and  eq  48  serves  as  a  definition  of  ice  strength  P*. 

The  two  main  unknowns  in  the  redistribution  theory  of  ridging  are  P(h) 
and  Y(hi,h2).  Thorndike  et  al.  (1975)  have  suggested  that  P(h)  be  given  by 

P(h)  =  max  [(l  -  g(h)dh/ci),0l  (49) 

where  ci  is  a  constant,  say  0.15.  The  basic  idea  here  is  that  under  ridg¬ 
ing  conditions  the  closure  of  both  thin  ice  and  thick  ice  occurs.  However, 
the  thinnest  ice  is  allowed  to  be  removed  to  a  greater  extent.  By  fixing 
ci  to  be  0.15,  only  the  thinnest  15%  of  the  ice  is  deformed.  Note  that  if 
there  is,  say,  25%  open  water,  only  open  water  will  be  removed  and  no  ridg¬ 
ing  will  occur. 

The  other  unknown  is  Y(hi,h.2),  which  specifies  how  ice  is  transferred 
from  one  category  to  another  by  ridging.  Thorndike  et  al.  (1975)  have  sug¬ 
gested  that 


Y(hi,h2)  *■  «(h2  -  Kh i ) K“ 1  (50) 
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where  K  is  a  constant.  Equation  50  effectively  states  that  ice  is  trans¬ 
ferred  into  a  fixed  multiple  of  its  own  thickness.  This  form  is  computa¬ 
tionally  simple  and  is  useful  for  examining  the  behavior  of  the  ice  thick¬ 
ness  distribution  model. 

However,  the  redistribution  of  eq  50  ignores  two  important  physical 
features  of  the  ridging.  One  feature  is  that  under  deformation,  ridging 
transfers  ice  to  a  variety  of  thickness  categories.  This,  for  example,  is 
evident  from  field  observations  of  ridges  (see,  e.g. ,  Weeks  et  al.  1971, 
Kovacs  1972)  ,  which  show  them  to  be  roughly  triangular.  These  observations 
indicate  that  under,  deformation,  leads  containing  thin  ice  of  a  given 
thickness  typically  form  ridges  having  a  triangular  cross  section.  To 
satisfy  eq  50  such  ridges  would  have  to  have  vertical  sides. 

The  second  physical  feature  is  that  typical  ridge  heights  divided  by 
the  thickness  of  ice  being  ridged  appear  to  decrease  with  increasing  thick¬ 
ness  (see,  e.g..  Tucker  and  Govoni  1981).  This  particular  scaling  is  an 
important  feature  of  the  Parmer  ter  and  Coon  (1972)  mechanistic  ridge  model. 
In  particular,  Parmerter  and  Coon  (1972)  found  that  simulated  ridges  using 
their  mechanistic  model  tended  to  have  a  maximum  limiting  height.  This 
limiting  height,  however,  tended  to  increase  more  slowly  than  the  thickness 
of  ice  being  ridged. 

To  avoid  some  of  these  shortcomings  Hibler  (1980a)  proposed  a  redis¬ 
tributor  which  distributed  ice  over  a  variety  of  thickness  categories  and, 
perhaps  more  importantly,  the  limiting  height  was  taken  to  scale  as  the 
square  root  of  the  thickness  of  ice  being  ridged.  Hibler  (1980a)  further 


noted  that  such  a  square  root  dependence  is  consistent  with  one's  intuition 
concerning  ridging.  Consider,  for  example,  two  equal  width  leads  under¬ 
going  ridging.  Assuming  that  all  the  thin  ice  is  deformed  into  piles  of 
relatively  small  blocks  with  similar  slope  angles,  the  two  ridges  formed 
would  have  heights  scaling  with  the  thickness  of  ice  being  ridged.  If  the 
distribution  of  lead  widths  is  independent  of  the  ice  thickness  in  leads, 
such  an  intuitive  argument  should  have  application. 

A  comparison  of  redistributors.  To  examine  the  ramification  of 
different  redistributors  let  us  compare  eq  50  to  several  other  forms  for 
Y(h i,h2)  that  have  been  suggested,  or  are  introduced  here  for  pedagogical 


purposes: 

Redistribution 

Description 

Y(h i,h2)  »  6(h2-Khi)  ^ 

(51) 

Thorndike  et  al.  (1975) 

with  K  *  constant 

Y(hi,h2) 


L 


2(H*-hi) 
with  H*  =  constant 


for  2hi<h1<2^  v'hT  (52)  Hibler  (1980a) 


Y(hj ,h2) 


2 

rr 


h  i(K  -4) 
with  K  *  constant 


for  2hi<h2<Kh i 


(53) 


Modified  Hibler  (1980a) 


Y(hi,h2)  *  6(h2-H') 


(54) 


Rubble  redistribution 


The  scaling  characteristics  of  different  redistributors  can  be  tested 
by  field  examination  of  ridge  height  and  blodt  size  characteristics.  A 
particularly  useful  data  set  for  this  purpose  has  been  obtained  and  analyz¬ 
ed  by  Tucker  and  Govoni  (1981).  This  data  set  was  taken  off  the  North  Slope 
of  Alaska.  Figure  13  shows  the  salient  characteristics  of  the  data  set 


Figure  13.  Compari¬ 
son  of  different  re¬ 
distribution  scaling 
laws  with  observed 
data.  The  solid 
curve  represents  the 
Hibler  (1980)  scal¬ 
ing  (eq  52) ,  the 
large  dashed  curve 
the  Thorndike  et  al. 
(197  5)  scaling  (eq 
51),  and  the  small 
dashed  curve  the 
"rubble”  scaling. 


V  /  V 


relevant  to  the  redistribution  theory.  The  error  bars  are  the  standard  er¬ 
ror  in  the  mean  ridge  height  estimate.  As  can  be  seen  from  this  figure  the 
square  root  scaling  fits  these  limited  data  better. 

To  approximately  convert  the  ridge  height  scaling  to  an  ice  thickness 
scaling,  a  4:1  ridge  keel  to  sail  ratio  is  assumed  (Kovacs  and  tfellor  1974). 
With  this  scaling  the  square  root  case  is  equivalent  to  an  H*  of  100  m  in 
the  Hibler  redistributor,  whereas  the  linear  scaling  represents  a  modified 
Hibler  redistributor.  The  rubble  case  represents  a  fixed  thickness  H*  -  20 
m. 

Of  particular  relevance  to  ice  mechanics  are  the  ice  strengths  gener¬ 
ated  by  different  redistributors.  To  get  some  feeling  for  these  strengths, 
consider  the  special  case  of  only  one  thickness  of  ice,  say  h0,  being 
ridged.  Formally  representing  this  condition  by  P(h)  -  6(h-ho)  one  can  sub¬ 
stitute  eq  46  in  eq  48  and  obtain  analytical  results  for  the  ice  strength. 
After  some  algebra  different  redistributors  yield  the  following  strength 
equations: 

ch2  (4 r3  -  3 r2  -i) 

^  3(  r-i>  r 

where  T  «  /H*/h0  for  the  Hibler  redistributor  and  T  *  K/2  for  the  modified 
Hibler  case, 

P*  -  Ch^K  (56) 

for  the  Thorndike  et  al.  (1975)  redistribution, 

P*  -  Ch0H'  (57) 

for  the  rubble  redistributor. 

3  3 

Using  a  representative  value  for  C  of  0.8026  x  10  N  m~  the  various 
scalings  in  Figure  13  yield  the  strength  versus  thickness  results  shown  in 
Figure  14.  In  the  constant  maximum  cutoff  case  we  have  used  the  modified 
Hibler  redistributor.  However,  from  the  above  strength  equation  it  is  clear 
that  the  Thorndike  et  al.  (1975)  definition  would  give  about  the  same 
strength  versus  thickness  variation,  but  with  all  strengths  50%  larger  than 
the  modified  Hibler  (1980a)  redistributor. 

The  main  feature  illustrated  by  Figure  14  is  that  as  expected  the 
square  root  scaling  tends  to  give  higher  strengths  for  thinner  ice  than  the 
constant  scaling.  The  reverse  is  true  for  thick  ice.  However,  in  both 
these  cases  the  three-dimensional  stress  will  increase  as  the  ice  becomes 
thicker.  This  is  in  contrast  to  the  rubble  case  where  the  three-dimensional 
stress  will  be  a  constant  at  ~1.6  *  104  N  m"2  (or  2.3  psi). 

Selected  seasonal  numerical  results.  The  Hibler  (1980a)  redistributor 
(eq  52)  has  been  used  as  one  component  of  a  coupled  variable  thickness  dy¬ 
namic  thermodynamic  sea  ice  model.  This  model  has  been  used  to  carry  out  a 
seasonal  equilibrium  simulation  of  the  Arctic  Basin  ice  cover  by  integrating 
the  model  for  five  years  with  the  same  seasonal  forcing  being  repeated  each 
year.  Of  particular  interest  to  this  review  are  the  ridge  buildup  results 
and  the  sensitivity  of  the  model  to  different  assumed  frictional  energy  dis- 
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Figure  14.  Ice  strength  versus 
thickness  for  different  redistribu¬ 
tors.  The  solid  and  long-dashed 
lines  are  for  a  uniform  redistribu¬ 
tion  with  the  maximum  thickness  scal¬ 
ing  as  the  square  root  (long-dashed 
line)  of  the  thickness  of  ice  being 
ridged  (eq  52),  and  linearly  (solid 
line)  with  the  thickness  of  ice  being 
ridged  (eq  51  and  53).  The  short- 
dashed  line  is  for  the  "rubble"  re¬ 
distributor  (eq  54)  where  all  ice  is 
transferred  into  a  20-m  thickness. 

In  these  curves,  frictional  losses 
are  assumed  to  be  equal  to  changes  in 
gravitation  potential  energy. 


sipation  values  in  the  ridging  process.  With  regard  to  ridge  buildup, 
while  the  model  does  not  explicitly  keep  track  of  ridging  it  is  possible  to 
determine  the  volume  of  deformed  ice  created  over  an  annual  cycle  at  a 
fixed  location.  Figure  15  shows  fifth-year  contours  of  this  quantity  using 
the  H*  value  of  50  m  in  eq  52  and  a  C  value  in  eq  48  of  2Cb  (see  eq  44). 

The  general  shape  of  the  roughness  contours  agrees  well  with  surface  rough¬ 
ness  observations  in  the  western  Arctic  Basin  (Weeks  et  al.  1971,  Hibler  et 
al.  1974).  The  two  major  features  are  a  heavy  buildup  of  ridging  off  the 
Canadian  Archipelago  and  less  ridging  in  the  Beaufort  Sea  than  near  the 
pole.  Note,  however,  that  there  is  a  zone  of  heavy  ridging  just  off  the 
North  Slope,  in  agreement  with  observations  by  Tucker  et  al.  (1979).  (This 
feature  is  masked  in  Figure  15  by  the  large  contouring  interval.  Finer 
contouring  intervals  show  a  peak  in  the  volume  of  deformed  ice  per  unit 
area  between  the  two  0.2-m  contours  in  Figure  15  off  the  North  Slope.)  The 
other  major  feature  is  a  tongue  of  high  ridging  further  offshore  near  the 
tip  of  Greenland.  Roughness  data  obtained  from  submarine  sonar  profiles  by 
LeShack  (private  communication)  also  show  such  a  tongue. 
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With  regard  to  ice  strengths  on  the  geophysical  scale,  a  major  unknown 
factor  is  the  amount  of  frictional  energy  loss  occurring  during  ridging. 

To  get  an  estimate  of  how  increases  in  this  assumed  frictional  loss  affect 
drift  simulations,  Hibler  (1980a)  also  carried  out  a  sensitivity  study  with 
the  frictional  energy  losses  assumed  to  be  9  times  the  potential  energy 
change  (C  «  10  C^) •  This  high  strength  case  resulted  in  an  unrealistic 
stoppage  of  flow  in  April  and  May  in  the  Arctic  Basin.  With  the  frictional 
loss  equal  to  the  potential  energy  change,  on  the  other  hand,  excessively 
large  ice  drifts  were  obtained.  Some  of  the  differences  between  these  two 
simulations  are  illustrated  in  Table  2  which  compares  average  net  simulated 
drift  of  three  ice  stations  obtained  during  a  335-day  period  of  the  simula¬ 
tion.  The  simulated  net  drift  values  were  obtained  by  interpolating  ice 
velocity  vectors  to  the  location  of  the  ice  stations  on  each  day  and  then 
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Figure  15.  Contours  (in  meters)  of  simulated  volume  of  deformed 
ice  per  unit  area  for  one  annual  cycle.  In  the  simulation  a  ten- 
level  ice  thickness  distribution  was  used  in  conjunction  with  a 
viscous-plastic  rheology.  A  description  of  the  model  is  given  by 
Hibler  (1980). 


Table  2.  Observed  and  simulated  ice  station  drift. 


Standard  simulation 

High  strength 

(C  =  2Cb) 

(C  =  10  Cb) 

Observed 

Net  ice  station  drift 

702  km 

339  km 

562  km 

Net  drift  angle 

26°  W 

13°  W 

6°  W 

multiplying  by  the  total  time.  In  addition  to  changing  C,  it  is  also  pos¬ 
sible  to  increase  the  strengths  by  increasing  H*  (see  eq  46).  With  respect 

to  maximum  P*  values  near  the  pole  in  these  two  s imula^ions, ^the  high 

strength  case  yielded  a  maximum  P*  value  of  abouj:  7*10^  N  m~  ,  whereas  the 

standard  case  yielded  a  maximum  P*  of  about  lx10  N  m"  .  Both  these  maxima 

occurred  in  March. 

Two-level  formulation  of  ice  strength 

The  multi-level  parameterization  provides  a  consistent  way  of  treating 
a  variable  thickness  ice  cover.  However,  as  is  apparent  from  the  above 
discussion  it  tends  to  transfer  much  of  the  complexity  of  the  problem  to 


the  mechanical  redistribution  function,  which  is  rather  complicated  and  in¬ 
herently  contains  a  number  of  assumptions.  The  idea  in  the  two-level  ap¬ 
proach  (Hibler  1979)  is  to  approximate  the  whole  thickness  distribution  by 
two  categories  of  ice:  thick  and  thin.  When  this  is  done  the  treatment  of 
the  ridging  process  becomes  much  simpler.  The  strength  is  then  taken  to  be 
an  empirical  function  of  these  two  categories,  with  the  most  important 
strength  dependence  being  on  the  thin  ice.  This  empirical  function  does 
contain  some  assumptions,  but  as  discussed  above  the  mechanical  redistribu¬ 
tion  function  also  contains  a  large  number  of  assumptions.  Also,  a  two- 
level  arctic  sea  ice  simulation  discussed  in  the  next  section  yields  re¬ 
sults  for  ice  thickness  and  drift  in  the  Arctic  Basin  very  similar  to  those 
obtained  from  the  multi-level  model,  which  indicates  that  for  many  climatic 
studies  the  two-level  approach  may  be  quite  useful. 

The  basic  concept  in  the  two-level  approach  is  that  many  of  the  dynam¬ 
ic  and  thermodynamic  characteristics  of  sea  ice  are  dominated  by  the  thin 
ice.  Consequently,  for  many  purposes  the  ice  thickness  distribution  may  be 
approximately  characterized  by  breaking  it  into  two  parts  (thick  and  thin). 
The  idea  here  is  that  if  the  thin  ice  component  is  monitored,  then  uncer¬ 
tainties  in  the  remaining  ice  thickness  distribution  may  not  be  critical. 
Within  this  two-level  approach  the  ice  cover  is  broken  down  into  an  area  A 
(often  called  the  compactness),  which  is  covered  by  thick  ice,  and  a  re¬ 
maining  area  1-A,  which  is  covered  by  thin  ice,  which,  for  computational 
convenience,  is  always  taken  to  be  of  zero  thickness  (i.e.  open  water). 

The  idea  here  is  to  have  the  open  water  approximately  represent  the  com¬ 
bined  fraction  of  both  open  water  and  thin  ice  up  to  some  cutoff  thickness 
ho.  The  remainder  of  the  ice  is  distributed  arbitrarily.  However,  since 
the  thin  ice  mass  is  normally  small,  the  mean  thickness  of  the  remaining 
"thick''  ice  is  approximately  equal  to  h/A. 

For  the  mean  thickness  h  and  compactness  A  the  following  continuity 
equations  are  used: 


3h/3t  »  -3(uh)/3x  -  3(vh)/3y  +  Sh  (58) 

3A/3t  >■  -3(uA)/3x  -  3(vA)/3y  +  (59) 

where  A  <  1,  u  is  the  x  component  of  the  ice  velocity  vector,  v  the  y  com¬ 
ponent  of  the  ice  velocity  vector,  and  and  are  thermodynamic  terms. 
While  eq  58  is  a  simple  continuity  equation  for  ice  mass  (characterized  by 
the  mean  thickness  h) ,  with  thermodynamic  source  and  sink  terms,  eq  59  is 
somewhat  more  complex.  By  including  the  restriction  that  A  <  1,  a  mechani¬ 
cal  sink  term  for  the  areal  fraction  of  ice  has  been  added  to  a  simple  con¬ 

tinuity  equation  for  the  ice  concentration.  This  sink  term  turns  on  when  A 
■  1  (i.e.,  no  open  water  left)  and  under  converging  conditions  removes 
enough  ice  area  through  ridging  to  prevent  further  increase  in  A.  Although 
the  sink  term  does  not  change  the  ice  mass,  it  can  cause  the  "thick"  ice 
thickness  to  increase  by  allowing  h  to  increase  while  A  does  not.  This 
simple  sink  term  is  the  equivalent  of  the  Thorndike  et  al.  (1975)  mechani¬ 
cal  redistribution.  Hence,  within  the  two-level  model  the  whole  ridging 
process  is  approximated  by  specifying  that  A  <  1.  Note  that  however  this 
ridging  problem  is  formally  treated  in  this  two-level  model  it  does  not 
affect  the  conservation  of  ice  mass,  which  is  explicitly  guaranteed  by  eq 
58.  This  explicit  treatment  of  the  ice  mass  continuity  equation  directly 
is  one  of  the  main  simplifications  of  the  two-level  model  over  the  multi¬ 
level  treatment. 
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Treatment  of  the  thermodynamic  terms  and  in  eq  58  and  59  varies 
with  investigators,  and  can  make  a  significant  difference  in  the  physics. 

In  ice  forecasting  applications  (e.g.  Nikiforov  1957,  Doronin  1970),  the 
term  is  often  omitted  and  A  is  simply  taken  as  the  area  covered  by  ice  (as 
opposed  to  open  water).  In  longer  term  seasonal  simulations  (Hibler  1979, 
Parkinson  and  Washington  1979)  has  been  used  to  cause  the  thin  ice  frac¬ 
tion  to  disappear  in  proportion  to  the  ice  growth  rates.  Under  these  for¬ 
mulations  (1-A)  represents  the  combined  fraction  of  thin  ice  and  open 
water. 

To  couple  the  ice  strength  to  the  ice  thickness  characteristics  the 
ice  strength  P*  is  taken  to  be  a  function  of  compactness  and  thickness  ac¬ 
cording  to 

P*  =  Pjjh  exp  [-C(l  -  A)]  (60) 

where  PQ  and  C  are  fixed  empirical  constants  and  h  is  in  meters.  This 
formulation  makes  the  strength  strongly  dependent  on  the  amount  of  thin  ice 
(characterized  by  the  compactness  A),  while  also  allowing  the  ice  to 
strengthen  as  it  becomes  thicker. 

While  directly  intuitive,  the  strength  expression  in  eq  60  can  also  be 
qualitatively  related  to  the  strength  assumptions  used  in  the  multi-level 
thickness  approach  discussed  above.  There,  under  converging  conditions, 
the  closure  of  both  thin  and  thick  ice  (and  also  open  water)  occurs  simul¬ 
taneously.  The  relative  amount  of  thin  ice  being  deformed,  however,  is 
weighted  by  the  area  of  thick  ice.  Hence,  for  small  amounts  of  thin  ice, 
more  thick  ice  is  deformed,  and  thus  the  ice  has  a  high  strength  which  is 
dependent  on  the  thick  ice  characteristics.  For  a  large  amount  of  thin 
i«_e,  on  the  other  hand,  mostly  thin  ice  is  deformed,  which  yields  low 
strength.  These  features  are  represented  in  eq  60  by  the  exponential 
dependence  on  (1  -A),  which  characterizes  the  amount  of  thin  ice  being 
ridged,  and  the  proportionality  to  h,  which  approximates  the  strength  of 
thick  ice  undergoing  deformation.  The  assumption  in  the  h  dependence  is 
that  the  distribution  of  thick  Ice  is  reasonably  well  described  by  its  mean 
(h/A).  However,  while  qualitatively  correct,  this  two-level  ice  thickness 
scheme  does  not  allow  the  particular  thickness  of  ice  involved  in  ridging 
to  be  determined. 

A  nunerical  comparison  of  thickness-strength  coupling 

To  give  some  feeling  for  the  role  of  thickness-strength  coupling  in 
the  two-level  and  multi-level  models  we  briefly  discuss  here  the  results  of 
simulations  of  the  Arctic  Basin  ice  cover.  In  order  to  compare  the  behav¬ 
ior  of  the  two-level  and  multi-level  models,  one-year  simulations  employing 
the  same  forcing  of  the  Arctic  Basin  ice  cover  are  compared.  These  simula¬ 
tions  consisted  of  integrating  the  coupled  continuity  and  momentum  balance 
equations  for  each  model  over  one  seasonal  cycle  for  the  Arctic  Basin  ice 
cover.  While  not  to  seasonal  equilibrium,  the  results  are  useful  for  com¬ 
paring  ice  strength  formulations.  Complete  details  on  the  numerical  frame¬ 
work  for  these  simulations  and  the  forcing  fields  are  described  in  Hibler 
(1979,  1980a).  Since  our  main  concern  is  the  ice  buildup  we  will  not  go 
into  the  details  of  the  forcing  fields  (and  heat  storage  in  the  oceanic 
boundary  layers),  which  are  discussed  in  the  above  references.  Th|  two 
main  parameter  sets  of  interest  for  this  discussion  are  PQ  (=  5*10  N  m“  )  , 


C  -  20,  the  strength  constant  for  the  two-level  case,  and  H*  (*50  m)  ,  the 
ridging  specification  parameter  for  the  multi-level  case  (see  eq  52).  In 
the  two-level  case  the  ice  strength  is  given  by  eq  60,  while  in  the  multi¬ 
level  case  the  strength  is  specified  by  combining  eq  48  with  eq  44,  46,  49 
and  52.  In  all  cases  a  uniform  initial  thickness  of  3.3  m  was  specified 
everywhere. 


Figure  16  shows  the  December  thickness  contours  for  both  models.  The 
boundaries  employed  in  the  simulation  are  specified  by  the  solid  lines. 
Also,  at  the  Greenland-Spitsbergen  passage  an  "open  boundary”  condition 
(see  Hibler  1979)  was  used,  which  allowed  ice  inflow  and  outflow.  As  can 
be  seen,  both  the  multi-level  and  two-level  models  yield  very  similar  ice 
thickness  buildup,  with  thicker  ice  near  the  archipelago  and  thinner  ice 
near  the  pole  and  off  the  Alaskan  and  Siberian  coasts.  This  thickness 
buildup  is  due  to  the  ice  dynamics,  which  tends  to  pile  up  ice  off  the 
archipelago  while  removing  it  from  other  regions.  In  both  models,  as  this 
ice  gets  thicker  it  gets  stronger  and  tends  to  prevent  further  convergence. 
The  main  feature  illustrated  by  this  figure  is,  however,  the  similarity  of 
the  thickness  contours  in  the  two  simulations,  despite  the  differences  in 
ice  thickness  parameter izations.  Both  these  thickness  contours  are  in  good 
agreement  with  observations  of  thickness  characteristics  in  the  Arctic  Ba¬ 
sin,  which  show  a  greater  buildup  along  the  archipelago,  with  typical  mean 
thicknesses  of  about  5-7  m  (see,  e.g.,  Wadhams  1981). 

It  should  be  noted  that  one  quantitative  difficulty  with  both  these 
simulations  is  the  use  of  eight-day  averaged  winds  (from  May  62  to  May  63) 
in  the  forcing  fields.  More  recent  two-level  simulations,  with  daily  time- 
varying  winds  of  the  Arctic  ice  cover  (Hibler  and  Walsh  1982)  and  portions 
of  the  Antarctic  ice  cover  (Hibler  and  Ackley  1982,  1983),  have  yielded 
higher  best  fit  strength  constants,  P0,  for  the  two-level  model  of  2.75*104 
N  m“  .  This  value  is  about  five  times  larger  than  the  best  fit  P*  value  use 
by  Hibler  (1979)  and  also  used  above.  The  main  reason  for  this  increase  is 
that  due  to  the  nonlinear  wind  stress,  daily  wind  fields  yield  a  larger 
average  wind  force.  As  a  consequence,  a  larger  ice  interaction  is  needed  to 
yield  modification  of  the  ice  drift  in  agreement  with  observations. 

Small  scale  ramifications  of  coupling  between 
ice  strength  and  thickness  characteristics 

As  discussed  above,  the  main  effect  of  thickness-strength  coupling  on 
the  large  scale  is  for  the  ice  to  become  stronger  in  regions  of  conver¬ 
gence,  which  in  turn  prevents  further  convergence.  While  this  same  behav¬ 
ior  also  occurs  on  smaller  scales,  there  can  be  a  number  of  fluctuating 
effects,  such  as  kinematic  waves,  which  are  masked  on  larger  scales.  Some 
of  these  effects  depend  on  the  nonlinear  ice  interaction  and  the  nonlinear 
coupling  between  ice  strength  and  thickness  characteristics.  To  yield  some 
Insights  into  these  small  scale  effects  Hibler  et  al.  (1983)  analyzed  some 
special  idealized  cases  of  the  two-level  model  using  the  same  grid  as  shown 
in  Figure  7.  The  essential  feature  of  these  results  is  that  for  a  nonlin¬ 
ear  model  it  is  possible  for  ice  velocities  at  a  given  location  to  fluctu¬ 
ate  in  time,  even  though  the  wind  forcing  may  be  fixed. 

Analytical  analysis  of  coupled  one-dimensional  systems.  To  demon¬ 
strate  the  overall  nature  of  plastic  buildup  on  a  small  scale  it  is  useful 
to  analyze  an  Idealized  one-dimensional  system.  As  in  the  moment um-bal- 


a.  Results  from  a  multi-level  simulation  employing  eight 
thickness  levels  and  an  ice  strength  based  on  an  assumed 
redistribution  process  (see  eq  44,  46,  48,  49  and  52). 
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b.  Results  from  a  two-thickness  level  simulation  using 
an  empirical  ice  strength  formulation  (eq  60). 

Figure  16.  Simulated  30-day  (day  331-360)  averaged  ice  thickness  contours 
in  meters  for  comparative  one-year  simulations. 
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Figure  17.  Geometry  of 
idealized  one-dimensional 
compactness  conditions. 


Figure  18.  Qualitative  plot  of 
velocity  versus  time  for  two  points 
in  idealized  one-dimensional  build¬ 
up  example. 


ance-only  case  (see  p.  11-14)  consider  a  one-dimensional  system  having  a 
linear  water  drag,  rigid  plastic  interaction,  and  constant  wind  field  t 
from  left  to  right.  To  simplify  analysis  the  strength  P*  is  taken  to 
depend  only  on  the  compactness  A.  The  initial  condition  consists  of  a  com¬ 
pactness  (A)  of  0.9.  For  the  coupling  of  P*  to  A  let  us  consider  P*  dis- 
continuously  changing  at  A  =  1  as  follows: 


P*(A) 


-Pi  if  A  <  1 
-P  if  A  *  1 


(61) 


where  P  >  Pj.  Because  of  this  discontinuous  change  the  local  momentum  bal¬ 
ance  will  not  be  modified  by  any  changes  in  compactness  until  the  compact¬ 
ness  goes  to  1.  Based  on  the  analysis  of  the  momentum  balance  only  in  sec¬ 
tion  4-c  all  deformation  takes  place  at  a  right-hand  boundary.  Consequent¬ 
ly,  at  some  arbitrary  time  At  after  the  initiation  of  motion  we  would  ex¬ 
pect  the  geometry  in  Figure  17  to  apply.  By  the  earlier  momentum  balance 
analysis,  for  (P-Pi)  <  H,3  all  velocities  in  the  L3  region  (compactness  of 
1)  will  be  zero.  Also  by  the  momentum  balance  analysis  the  velocities  in 
the  L2  region  will  be  given  by 

U2  -  (tL2  -  Pi)/cL2  (62) 

so  that  as  L2  becomes  smaller  U2  will  slow  down  somewhat  unless  Pj  ■  0. 

Since  with  the  L3  region  motionless  all  the  deformation  is  occurring  at  the 
L2,L3  boundary,  it  is  clear  that  the  L2,L3  boundary  will  propagate  outward 
(to  the  left)  at  a  speed  of  U2/0.1.  However,  note  that  for  L 3  large  enough 
(i.e.  L3  >  ( P-P 1 ) / t)  the  L3  region  will  start  moving  again  with  a  speed 

U3  =  (TL3  “  (P-Pi))/cL3.  (63) 

When  this  occurs  the  propagation  of  the  ridging  front  will  slow  down  some¬ 
what  to  (U2-U3)/0.1.  Also,  once  all  the  ice  is  converted  to  100%  compact¬ 
ness,  the  speed  will  then  slowly  decrease  as  the  length  L3  decreases  due  to 
removal  of  ice  by  ridging  at  the  right-hand  boundary.  Based  on  these  con¬ 
siderations  one  can  construct  a  qualitative  speed  versus  time  of  two  arbi¬ 
trary  points,  say  p  and  p’ ,  where  p’  is  located  further  left  than  p  (Fig. 
18). 


Obviously,  by  varying  strengths  and/or  the  wind  forcing  the  various 
slopes,  magnitudes  and  spacing  in  time  of  these  profiles  can  be  modified. 
One  of  the  more  interesting  features  of  the  system  is  that  due  to  the  com¬ 
plexities  of  the  thickness  coupling,  for  appropriate  strength  differences, 
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points  near  the  right-hand  coast  will  actually  stop  for  a  while  and  then 
speed  up,  even  though  external  forcing  is  constant.  The  essential  physics 
here  is  that  because  of  the  discontinuous  nature  of  the  ridging  front  and 
the  nonlinear  rheology,  points  in  the  L3  region  will  not  move  at  all  until 
L3  becomes  large  enough  to  accumulate  an  adequate  wind  fetch  to  overcome 
the  differential  of  plastic  stresses  at  the  boundaries. 

Numerical  Investigation  of  response  characteristics  of  coupled  system 
A  series  of  two-dimensional  numerical  simulations  yielded  results  that 
generally  agree  with  the  one-dimensional  analysis.  However,  there  are 
significant  differences  associated  with  the  fact  that  P*  does  not  change 
discon tinuously  but  is  instead  a  smooth  function  of  A.  For  the  standard 
case  study,  the  coupled  equations  were  integrated  for  18  hours  at  15-minute 
time  steps  using  the  grid  of  Figure  7.  The  initial  conditions  consisted  of 
10%  open  water  together  with  a  mean  ice  thickness  of  0.5  m.  To  simplify 
analysis  the  ice  strength  is  taken  to  depend  only  on  the  compactness,  not 
the  ice  thickness.  As  a  consequence  the  ice  strength  in  these  particular 
idealized  studies  P*  is  given  by 

P*  -  (0.5  ■)  P  e-C(1-A)  (64) 

o 

To  make  the  test  more  comparable  to  the  one-dimensional  analysis,  turning 
angles  and  the  Coriolis  parameter  were  set  equal  to  zero.  A  constant  wind 
speed  of  9.23  m/s  in  the  positive  x  direction  was  used,  together  with  an 
ice  strength  constant  PQ  of  2. 0*10 4  N  m"  ,  and  a  constant  C  of  30.  Other 
parame  ters  are  the  same  as  in  Table  1 . 

The  basic  characteristics  of  this  idealized  numerical  simulation  are 
shown  in  Figure  19,  which  shows  x  velocity  components  versus  time  for  sev¬ 
eral  grid  points.  These  curves  were  taken  from  grid  points  and  grid  cells 
centered  in  the  y  direction.  The  general  behavior  of  the  velocity-time 
series  is  commensurate  with  the  idealized  one-dimensional  analysis,  with 
the  points  nearest  the  coast  rapidly  decreasing  in  speed,  with  a  later 
increase  as  the  region  of  ridging  moves  further  out.  However,  there  are 
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Figure  19.  Time  series  of  x  component  of  ice  velocity 
progressively  further  from  right-hand  boundary  of  numer¬ 
ical  buildup  experiment.  The  grid  points  plotted  were 
centered  in  the  y  direction  and  the  distances  from  the 
right-hand  boundary  are  labeled. 


significant  differences  from  the  analytic  case  associated  with  the  fact 
that  the  motion  near  the  coast  is  balanced  by  a  gradient  in  compactness 
rather  than  a  discontinuous  change.  This  gradient  means  a)  that  the  points 
do  not  ever  totally  stop  (since  some  convergence  oust  occur  everywhere  to 
allow  the  gradient  to  move  outward)  and  b)  that  the  kinematic  waves  are 
allowed  to  propagate  upstream  (or  more  precisely  down  the  stress  gradient). 

These  waves  arise  from  the  fact  that  the  velocity  of  ice  is  balanced 
in  part  by  the  stress  gradient,  while  the  stress  gradient  is  maintained  by 
a  convergence  of  the  ice  velocity  field.  As  a  consequence,  a  perturbation 
of  either  one  of  these  quantities  can  cause  a  wave  to  propagate  down  the 
stress  gradient  slope.  This  wave  can  naturally  develop  in  the  evolution  of 
the  Systran  since  the  initial  buildup  can  be  out  of  equilibrium  with  regard 
to  kinematic  waves  (when  the  dependence  of  strength  on  compactness  is  non¬ 
linear).  Such  a  kinematic  wave  is  apparent  in  Figure  19  at  hour  6  in  the 
velocity  time  series  of  the  point  nearest  the  coast.  Tests  at  smaller  time 
steps  verified  that  this  wave  was  not  a  numerical  artifact.  Closer  exam¬ 
inations  of  this  wave  show  it  to  slow  down  somewhat  and  diffuse  as  it  pro¬ 
gresses  outward.  More  detailed  investigations  of  the  characteristics  of 
these  waves  are  needed. 


THE  ROLE  OF  ICE  DYNAMICS  IN  THE  ATMOSPHERE-ICE-OCEAN  SYSTEM 
General  discussion 

The  fact  that  ice  is  moving  and  deforming  rather  than  static  has  im¬ 
portant  ramifications.  A  feature  that  has  been  emphasized  a  great  deal  re¬ 
cently  is  the  fact  that  dynamics  creates  thin  ice  and  leads  in  an  ice  cover 
(Maykut  1978,  Washington  et  al.  1976).  However,  in  addition  to  creating 
leads  the  dynamics  also  causes  advectlve  effects  which  can  play  a  pivotal 
role  in  the  atmospheric  and  oceanic  circulation.  In  the  presence  of  a  free 
ice  edge,  for  example,  advectlve  effects  can  transfer  ice  to  the  marginal 
regions  to  be  rapidly  melted.  When  viewed  from  a  satellite  photo  such  an 
ice  cover  may  appear  static,  since  the  ice  extent  does  not  rapidly  change. 
However,  the  actual  physics  maintaining  such  a  ’’static"  extent  can  be  com¬ 
plex  and  largely  dependent  upon  processes  occurring  at  the  ice  margin. 

Advection  also  causes  net  imbalances  of  ice  growth  by  transferring  ice 
from  regions  of  high  growth  to  regions  of  greater  melt.  Since  most  of  the 
salt  is  taken  out  of  the  ice  when  it  freezes,  the  high  growth  regions  will 
be  regions  of  high  salt  flux  into  the  ocean,  whereas  the  high  melt  regions 
will  be  areas  of  fresh  water  input.  These  features  can  significantly  alter 
the  oceanic  circulation  and  vertical  convection.  With  regard  to  the  atmo¬ 
sphere,  regions  of  high  growth  are  regions  of  high  heat  exchange  with  the 
atmosphere  due  to  the  latent  heat  of  freezing  of  the  ice  cover.  A  third, 
rather  indirect  way  ice  dynamics  can  affect  the  atmosphere-ocean  heat  ex¬ 
changes  has  recently  been  noted  by  Hibler  and  Ackley  (1982).  They  carried 
out  several  simulations  invoking  an  ice  feedback  mechanism  similar  to  that 
used  in  many  simple  climate  models,  whereby  the  presence  of  ice  caused  the 
lower  atmosphere  to  be  cooled.  With  this  feedback  mechanism,  simulations 
without  dynamics  produced  very  unrealistic  summer  ice  margins,  whereas 
simulations  with  dynamics  produced  a  realistic  delay.  Basically,  the 
dynamics  affords  a  mechanism  for  "turning  off"  the  feedback  link,  which  is 
not  possible  in  the  purely  static  thermodynamic  models. 


Finally,  in  addition  to  modifying  the  thermal  and  salinity  fluxes  at 
the  air/ocean  interface,  sea  ice  can  substantially  modify  the  momentun 
transfer  by  wind  to  the  ocean.  Specifically,  if  we  consider  the  ice-ocean 
system  to  be  coupled,  then  the  momentum  transfer  into  the  ocean  is  given  by 
the  wind  stress  less  the  force  due  to  ice  interaction.  As  discussed  ear¬ 
lier,  due  to  the  highly  nonlinear  ice  interaction  the  ice  stress  terms  can 
be  very  complex. 

To  illustrate  the  relative  importance  of  some  of  these  features,  re¬ 
sults  from  two  series  of  simulations  of  the  Weddell  Sea  pack  ice  are  dis¬ 
cussed  in  the  following  subsections.  The  first  subsection  describes  simu¬ 
lations  of  the  Weddell  Sea  pack  ice,  while  the  second  subsection  examines 
the  salient  characteristics  of  an  ice-covered  Arctic  Ocean  model.  A  more 
complete  description  of  these  simulations  is  given  elsewhere  (Hibler  and 
Ackley  1983,  Hibler  and  Bryan  1984). 

Resul ts  from  a  Weddell  Sea  pack  ice  model 


Model  description.  The  model  essentially  consists  of  a  momentum  bal¬ 
ance  employing  a  plastic  rheology  coupled  to  the  two-thickness  evolution 
equations  for  the  ice  thickness  and  compactness  discussed  earlier  (eq  58- 
60).  These  equations  also  include  ice  growth  and  decay  due  to  surface  heat 
budget  considerations  and  heat  storage  in  a  motionless  fixed-depth  oceanic 
boundary  layer.  For  the  ice  strength  (see  eq  60)  we  take  PQ  -  2.75  x  10**  N 
m"  and  C  =  20.  This  PQ  value  is  about  five  times  larger  than  the  PQ  value 
used  by  Hibler  (1979).  This  increase  was  necessitated  by  the  fact  that 
daily  wind  fields  were  used  here  as  opposed  to  the  eight-day  average  wind 
fields  used  by  Hibler  (1979). 


Following  formulations  developed  by  Hibler  (1979),  the  thermodynamic 
terms  in  the  continuity  equations  (58  and  59)  are  given  by 


Sh  =  f (h/ A) A  +  (1  -  A)f(0) 


SA  = 


(f(0)/ho)(l  -  A) 

0 


if  f(0)  >  0 
if  f(0)  <  0 


(65) 


+ 


0 

(A/2h)Sh 


if  Sh  >  0 
if  Sh  <  0 


(66) 


with  f(h)  the  growth  rate  of  ice  of  thickness  h,  and  hQ  a  fixed  demarcation 
thickness  between  thin  and  thick  ice. 


The  Sf,  term  specifies  the  net  ice  growth  or  melt.  Within  the  two- 
level  approximation,  Sy,  jls  given  by  the  sum  of  the  ice  grown  on  open  water 
plus  the  additional  growth  over  the  portion  of  the  cell  covered  by  thick 
ice.  To  approximate  the  growth  and  decay  rate  of  this  thick  ice,  its  mean 
growth  rate  is  estimated  to  be  that  of  ice  of  constant  thickness  h/A.  For 
melting  conditions  the  same  sum  over  open  water  and  thick  ice  is  used.  A 
critical  assumption  here  is  that  the  heat  absorbed  by  open  water  will  hori¬ 
zontally  mix  and  melt  additional  ice  until  the  mixed  layer  returns  to 
freezing. 
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The  Sa  term  characterizes  the  way  in  which  growth  and  decay  change  the 
relative  areal  extents  of  thin  and  thick  ice.  The  basic  physical  notion 
embodied  by  this  term  is  that  the  areal  fraction  of  thin  ice  will  decrease 
rapidly  under  freezing  conditions,  and  increase  slowly  under  melting  condi¬ 
tions.  To  parameterize  the  freezing  effect  simply,  the  fraction  of  open 
water  (1  -  A)  is  allowed  to  decay  exponentially  with  a  time  constant  of 
ho/f(0),  which  gives  the  first  term  in  eq  66.  The  constant  h0  is  chosen  to 
be  small  compared  to  mean  ice  thicknesses  but  large  enough  so  that  heat 
fluxes  through  hQ-thick  ice  are  substantially  less  than  through  open  water. 
In  the  Arctic  simulations  performed  by  Hibler  (1979)  ho  ■  0.5  m  is  used. 
However,  for  these  Antarctic  simulations  ho  is  set  equal  to  1.0.  This  lar¬ 
ger  value  is  based  on  field  observations  reported  in  Ackley  et  al.  (1980) 
and  Gow  et  al.  (1982),  which  suggest  that  frazil  ice  formation  may  well 
prolong  ice  production  in  open  water  regions.  Other  than  ho  and  Pot  the 
numerical  parameters  used  in  the  simulation  are  identical  to  those  used  by 
Hibler  (1979)  and  Hibler  and  Walsh  (1982). 

The  second  term  in  eq  66  accounts  for  melting.  Its  magnitude  is  de¬ 
rived  by  assuming  that  the  thick  ice  is  uniformly  distributed  between  0  and 
2  h/A  in  thickness,  and  all  melts  at  the  same  rate.  Therefore,  over  a  time 
At  the  ice  of  thickness  less  than  At  will  melt  and  form  open  water.  By 
the  assumption  of  uniform  distribution  this  ice  covers  a  fraction  of  area 
equal  to  At  A/2h,  which  for  At  small  yields  the  second  term  in  the  equa¬ 
tion  for  Sa* 

The  vertical  growth  rates  (f(h))  are  estimated  by  employing  Semtnc*'s 
(1976)  sea  ice  thermodynamic  model  in  conjunction  with  a  surface  heat  bud¬ 
get  calculation  similar  to  that  of  Parkinson  and  Washington  (1979)  and 
Manabe  et  al.  (1979).  The  basic  surface  energy  budget  calculation  (for  de¬ 
tails  see  Hibler,  1980a,  Appendix  B)  Includes  incoming  shortwave  and  long¬ 
wave  radiation  and  sensible  and  latent  heat  fluxes,  and  provides  for  the 
determination  of  the  surface  temperature  of  the  ice  by  numerical  iteration. 

Following  Manabe  et  al.  (1979),  the  effects  of  snow  cover  are  approxi¬ 
mated  by  allowing  the  ice  surface  albedo  to  be  that  of  snow  when  the  calcu¬ 
lated  surface  temperature  is  below  freezing,  and  that  of  snow-free  ice  when 
the  surface  temperature  of  the  ice  is  at  the  melting  point.  The  model  also 
includes  a  motionless  oceanic  mixed  layer  of  fixed  depth  (30  m)  and  assumes 
a  constant  oceanic  heat  flux  of  2  W/m  .  In  the  two-level  simulation,  when 
open  water  is  absorbing  heat  all  of  the  absorbed  heat  is  used  to  further 
melt  the  ice.  If  there  is  no  ice  the  absorbed  heat  raises  the  temperature 
of  the  mixed  layer.  Under  growth  conditions,  no  ice  is  allowed  to  form 
until  the  mixed  layer  reaches  the  freezing  temperature  of  seawater. 

Input  fields.  Input  fields  to  drive  the  model  consist  of  average  dai¬ 
ly  surface  atmospheric  air  temperatures,  mixing  ratios,  and  surface  pres¬ 
sure  taken  from  the  1979  Australian  analysis.  Preliminary  simulations  with 
these  analyzed  atmospheric  temperature  fields  yielded  unrealistically  small 
ice  concentrations  in  summer.  Comparison  of  analyzed  surface  temperature 
with  buoy  observations  indicated  that  the  main  problem  was  the  excessively 
high  temperatures  in  the  late  summer  and  fall  (February-April)  in  the  Aus¬ 
tralian  surface  analysis,  which  did  not  Include  these  buoy  observations. 

Two  procedures  were  examined  to  construct  a  more  realistic  surface  analy¬ 
sis.  In  a  crude  reanalysis  method,  hereafter  referred  to  as  the  standard 
case,  air  temperatures  and  mixing  ratios  were  not  allowed  to  rise  above 
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Figure  20.  Grid  used  for  numerical  simulations  of  the 
Weddell  Sea  pack  ice.  The  dashed  line  encloses  the 
area  where  the  temperature  and  humidity  fields  were  re¬ 
analyzed.  The  hatched  grid  cells  represent  open  bound¬ 
aries  which  allow  inflow  or  outflow  of  ice,  and  are 
taken  to  have  0  strength  for  ice  dynamics  calculations 
(see  Hibler,  1979,  for  a  more  complete  discussion  of 
these  open  boundaries).  The  open  circles  denote  the  1 
January  1979  positions  of  two  drifting  buoys. 

271.2  K  and  0.25  g  kg-1  in  the  region  designated  by  the  square  box  in  Fig¬ 
ure  20  in  all  months  except  January.  Slightly  different  regions  were  ex¬ 
perimented  with  but  had  little  effect  on  the  results.  For  the  more  fully 
coupled  case,  the  air  temperature  and  the  humidity  were  related  to  the  pre¬ 
sence  of  ice  by  specifying  maxima  of  271.2  K  and  0.25  g  kg-1  for  the  tem¬ 
perature  and  mixing  ratios  if  the  ice  cover  concentration  was  greater  than 
0.50  (or  in  the  case  of  the  thermodynamics-only  model  if  the  ice  thickness 
was  less  than  25  cm).  It  is  notable  that  such  a  "feedback"  mechanism  is 
consistent  with  field  observations  (S.F.  Ackley,  private  communication) 
near  the  ice  edge,  suggesting  that  the  air  temperature  undergoes  an  abrupt 


change  as  the  ice  edge  is  crossed.  Also,  measurements  of  atmospheric 
boundary  layer  modification  (Andreas  et  al.  1979)  support  the  assertion  of 
a  feedback  between  atmospheric  temperature  and  ice  extent. 


For  the  calculation  of  geostrophic  ocean  current  fields,  mean  dynamic 
topography  values  reported  by  Gordon  et  al.  (1978)  were  used.  In  the  radi¬ 
ation  calculations,  >ar ame ter izat ions  similar  to  those  employed  by  Parkin¬ 
son  and  Washington  (1979)  were  used.  Specifically,  daily  global  solar  ra¬ 
diation  under  cloudless  skies  was  obtained  by  integrating  an  empirical 
equation  by  Zillmann  (1972)  over  solar  zenith  angles  for  any  particular 
day.  [Zenith  angles  at  half-hour  intervals  for  this  purpose  were  obtained 
from  a  numerical  solution  of  Kepler's  equation  (Holloway,  personal  communi¬ 
cation).]  Incoming  longwave  radiation  was  obtained  using  Idso  and  Jack¬ 
son's  (1969)  formula  for  clear  skies.  For  cloud  cover,  a  constant  value  of 
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0.85  was  assumed  at  all  points.  This  value  is  commensurate  with  maximum 
estimates  by  Van  Loon  (1972). 


Simulation  results.  Using  the  modified  input  fields  and  feedback  pa¬ 
rameterization  discussed  above,  two-year  simulations  were  carried  out  using 
a)  the  full  coupled  model  and,  for  comparison,  b)  only  the  thermodynamic 
portion  of  the  model,  and  c)  a  version  including  leads  but  no  advection. 
This  latter  case  was  obtained  by  using  the  continuity  equations 

•44  “  ( V  *u)A  +  S  (67) 

o  t  —  a 

-  Sh  -  f(h/A)A  +  (l-A)f(O).  (68) 


To  generate  the  leads,  the  ice  velocity  field  from  the  second  year  of  the 
standard  fully  coupled  simulation  was  used.  Note  that  no  ridging  is  con¬ 
sidered  in  this  case.  Basically,  with  this  "leads-only"  case  deformation 
only  enters  by  affecting  A,  which  in  turn  modifies  the  average  growth  rate 
Sh  over  a  given  region.  In  all  cases  the  initial  condition  consisted  of  2 
m  of  ice  everywhere.  Comparison  of  these  results  gives  some  insight  into 
the  role  of  ice  dynamics  on  the  ice  edge  fluctuations  and  the  effects  of 
ice-atmosphere  feedback  on  the  seasonal  cycle  of  ice. 

Some  selected  ice  extent  results  from  these  simulations  are  illustrat¬ 
ed  in  Figure  21,  which  shows  plots  of  the  ice  concentration  for  the  full 
model  and  thin  ice  contours  for  the  thermodynamics-only  case  employing  the 
standard  forcing.  The  observations  are  based  on  Fleet  Weather  ice  charts 
(prepared  by  the  Navy-NOAA  Joint  Ice  Center,  Suitland,  Md.)  and  denote  50% 
ice  concentration  limits.  In  most  cases  this  50%  limit  is  effectively  at 
the  ice  edge  estimated  in  the  Fleet  Weather  charts.  In  this  figure  two 
different  days  were  chosen  to  represent  salient  features  of  both  the  simu¬ 
lated  and  observed  ice  extents.  Day  18  (18  January)  represents  a  period 
close  to  the  end  of  the  rapid  decay  and  shows  a  characteristic  geometric 
tongue  of  ice  often  seen  at  this  time  of  year.  Day  298  (25  October),  on 
the  other  hand,  represents  a  typical  maximum  extent  at  a  time  just  prior  to 
the  beginning  of  the  rapid  decay  phase. 

As  can  be  seen  from  Figure  21,  on  day  18  the  full  model  reproduces  the 
tongue  effect  and  ice  extent  much  better  than  the  thermodynamics-only  case. 
Analysis  of  the  ice  velocity  field  (see  the  next  section)  shows  that  this 
tongue  is  largely  due  to  northward  and  eastward  advection  of  ice  by  the  ice 
velocity  field.  (The  tongue  is,  however,  a  bi  t  too  far  south,  partially 
due  to  insufficient  northward  ice  transport  as  discussed  in  the  next  sec¬ 
tion.)  In  the  thermo  dynamics -only  case  the  structure  of  the  ice  edge  is 
totally  incorrect  and  there  is  too  much  ice. 

On  day  298  the  results  are,  however,  quite  similar,  with  both  the  full 
model  and  the  thermodynamic  case  showing  extents  slightly  less  than  ob¬ 
served.  This  in  turn  indicates  that  the  expansion  to  maximum  is  largely 
thermodynamic  in  nature  and  controlled  (in  both  simulations)  by  the  temper¬ 
ature  field. 

The  effect  of  dynamics  on  the  seasonal  variation  is  illustrated  more 
graphically  in  Figure  22,  which  shows  the  observed  time  series  of  the  total 
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b.  Results  using  the  feedback  parameterization  between  the 
atmosphere  and  the  ice  discussed  in  the  text. 


Figure  22.  Time  series  of  total  area  covered  by  ice  within 
the  simulation  region. 
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grid  area  covered  by  ice  with  concentrations  greater  than  15%.  This  figure 
includes  simulations  using  both  the  "feedback"  (part  b)  forcing  and  the 
"standard”  forcing  (part  a).  Also  shown  in  this  figure  are  the  time  series 
of  the  thermodynamics-only  case  and  of  the  "leads-only"  simulation.  In  the 
thermodynamics -only  case  the  area  given  is  for  ice  of  non-zero  thickness. 

As  can  be  seen,  the  coupled  model,  with  the  dynamics  included,  gives  a  big¬ 
ger  seasonal  swing  of  extent  and  a  more  rapid  decay  of  the  ice  cover.  Both 
these  effects  are  primarily  due  to  the  large  amount  of  deformation  and  ad- 
vection  caused  by  the  dynamics,  which  in  turn  creates  open  wat^r.  As  the 
shortwave  radiation  terms  increase  during  spring  and  summer  th-i-S  open  water 
absorbs  large  amounts  of  heat  and  causes  a  rapid  decay  of  the  ice  cover. 
However,  a  difficulty  with  the  standard  case  (and  the  thermodynamics-only 
case)  is  a  lack  of  persistence  of  the  ice  cover  in  the  spring  (day  250  on). 

As  can  be  seen,  the  full  model  reproduces  a  realistic  seasonal  varia¬ 
tion  and  geographical  extent.  This  feature  is  substantially  improved  in 
the  feedback  case,  which  produces  the  persistence  in  the  ice  extent  in  the 
spring  (day  250  on)  that  is  in  slightly  better  agreement  with  observation 
than  the  standard  case  discussed  above,  without  significantly  modifying  the 
rapid  decay  seen  in  the  summer  (December-January)  period.  Although  not 
shown,  the  feedback  also  produces  a  realistic  ice  tongue  similar  to  the 
standard  case  discussed  previously.  The  thermodynamics-only  case,  on  the 
other  hand,  effectively  produces  a  constant  ice  edge  with  very  little  sea¬ 
sonal  cycle. 

With  respect  to  the  simulation  including  leads,  its  behavior  is  simi¬ 
lar  to  the  thermodynamics-only  case  except  that  it  has  a  bigger  seasonal 
swing.  However,  the  inclusion  of  leads  alone  without  advection  is  not  ade¬ 
quate  to  give  a  realistic  seasonal  swing.  Also,  examination  of  the  geo¬ 
graphical  contours  of  the  leads-only  case  shows  them  to  be  very  similar  to 
the  thermodynamics-only  case,  and  especially  unrealistic  in  summer.  Analy¬ 
sis  of  the  results  shows  that  in  the  full  model  the  reduced  extent  is  par¬ 
tially  due  to  winds  near  the  ice  edge  —  winds  which  tend  to  shift  the 
whole  ice  edge  eastward  and  thus  reduce  the  extent. 

A  notable  feature  of  these  simulations  is  that  for  the  forcing  em¬ 
ployed  here  it  is  not  necessary  to  increase  the  oceanic  heat  flux  above 
typical  Arctic  values  to  obtain  a  realistic  seasonal  cycle  when  realistic 
dynamics  are  included.  This  is  in  contrast  to  calculations  based  primarily 
on  thermodynamic  considerations  by  previous  authors,  which  have  suggested 
the  need  for  a  larger  oceanic  heat  flux  to  explain  the  ice  decay. 

While  there  are  differences  in  extent,  a  much  more  major  difference 
between  the  full  model  and  the  deformation-only  and  thermodynamics-only 
models  is  in  the  net  seasonal  ice  growth.  To  examine  these  characteristics 
the  fully  coupled  feedback  simulation  was  extended  to  a  third  year,  and  the 
net  annual  ice  growth  calculated.  Figure  23  shows  the  contours  of  this  net 
ice  growth  over  an  annual  cycle.  In  both  the  leads-only  and  thermodynam¬ 
ics-only  cases  the  net  ice  growth  is  very  small,  since  in  seasonal  equilib¬ 
rium  all  ice  growth  is  balanced  by  melt.  However,  with  advection  included 
there  are  drastic  changes,  with  much  of  the  ice  growth  taking  place  in 
southern  parts  of  the  grid  and  ice  melt  occurring  farther  north  near  the 
Equator.  Off  the  Fllchner  Ice  Shelf  the  ice  growth  exceeds  3  m.  High 
growth  in  this  region  is  known  to  be  a  major  factor  in  bottom  water  forma¬ 
tion,  and  the  numbers  obtained  from  the  simulation  are  commensurate  with 
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Figure  23.  Contours  (In  meters)  of  the  net  annual  ice  growth 
over  an  annual  cycle  for  the  fully  coupled  dynamic  thermo¬ 
dynamic  model.  The  contours  were  compiled  from  the  third  year 
of  the  feedback  simulation. 


maximum  estimates  in  the  literature  (see,  e.g. ,  Ackley  1980)  based  on 
oceanographic  considerations.  Also,  the  fact  that  ice  melting  takes  place 
at  different  locations  than  growth  means  that  large  amounts  of  fresh  water 
can  effectively  be  transported  to  different  regions  of  the  Weddell  Sea,  a 
factor  of  some  importance  in  the  pycnocline  structure  there.  Note  also 
that  much  of  the  melt  is  highly  concentrated  at  the  ice  edge,  illustrating 
a  conveyor-belt-like  effect  whereby  ice  is  transported  to  the  edge  to  be 
melted. 


Ice  edge  characteristics  of  an  ice-covered  Arctic  Ocean  model 


In  the  simulations  of  Antarctic  sea  ice  discussed  above,  realistic  ice 
edge  results  were  obtained  with  the  ocean  parameterized  by  only  a  motion¬ 
less  fixed-depth  mixed  layer.  In  the  Greenland/Norwegian  Seas,  however, 
such  an  approximation  (see,  e.g.,  Hibler  and  Walsh  1982)  yields  unrealisti¬ 
cally  large  ice  extents.  This  is  not  surprising  in  that  large  scale  ocean¬ 
ic  circulation  effects  would  be  expected  to  be  particularly  pronounced  in 
the  Greenland  and  Norwegian  Seas  due  to  warm  northward  currents  encounter¬ 
ing  rapidly  cooling  atmospheric  conditions  and  southward-advancing  sea  ice. 
In  order  to  examine  the  dominant  effects  of  the  full  three-dimensional 
ocean  circulation  on  the  simulated  ice  margin  and  to  examine  ice-ocean 
coupling  effects,  Hibler  and  Bryan  (1984)  constructed  a  diagnostic  ice- 
ocean  model  and  used  it  to  carry  out  a  series  of  seasonal  simulations  of 
the  Arctic,  Greenland  and  Norwegian  Seas.  The  term  "diagnostic”  refers  to 
the  fact  that  temperature  and  salinity  fields  are  largely  specified  by  ob¬ 
served  data  by  a  relaxation  method  rather  than  fully  simulated.  However, 
as  discussed  below  this  relaxation  is  done  in  a  relatively  weak  manner  so 
that  seasonal  and  month-to-month  variations  in  the  temperature  and  salinity 
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due  to  ice-ocean  coupling  effects  will  be  simulated  by  the  model.  Some  of 
the  salient  results  of  these  simulations  are  briefly  described  here. 

In  this  study  an  existing  dynamic- thermo  dynamic  sea-ice  model  (Hibler 
1979,  1980)  was  coupled  with  a  multi-level  baroclinic  ocean  model  (Bryan 
1969).  The  sea-ice  model  (essentially  equivalent  to  the  Weddell  Sea  pack 
ice  model  described  above)  supplied  heat  flux,  salt  flux,  and  momentum 
exchange  boundary  conditions  for  the  top  of  the  ocean.  The  ocean  model,  in 
turn,  supplied  current  and  heat  exchange  information  to  the  ice  model. 

Since  the  concern  here  was  to  examine  the  effect  of  ocean  circulation  on 
sea  ice,  observed  oceanic  temperature  and  salinity  data  (Levitus  1982)  were 
used  to  weakly  force  the  ocean  model  so  that  its  equilibrium  time  scale  is 
similar  to  that  of  the  ice  model  (3  to  5  years).  This  forcing  is  done  by 
adding  additional  source  terms  to  the  heat  and  salt  conservation  equations 
of  the  form  R  [T0  -  Tl  and  R  [S0  -  S].  Here  R  is  a  relaxation  constant 
taken  to  be  0.333  yr“  ,  and  T0  and  S0  are  the  mean  annual  observed  tempera¬ 
ture  and  salinity  fields.  This  "diagnostic"  method  allows  the  ocean  model 
to  be  forced  to  available  climatological  ocean  data,  while  at  the  same  time 
allowing  considerable  adjustment  in  the  upper  ocean  due  to  the  effects  of 
ice/ocean  interaction.  In  addition,  the  barotropic  mode  of  the  ocean  is 
fully  simulated  so  that  temporally  varying  currents  due  to  surface  stress 
fluctuations  are  part  of  the  model  predictions.  The  basic  idea  here  is  to 
specify  rather  than  model  the  ocean  circulation  over  climatic  time  scales 
while  modeling  the  shorter  term  fluctuations. 

The  most  noticeable  effect  of  the  modeled  ocean  circulation  is  to 
greatly  improve  the  ice  margin  simulation,  as  shown  in  Figure  24a.  This 
Improved  ice  margin  is  due  to  the  large  oceanic  heat  flux  from  the  deeper 
ocean  into  the  upper  mixed  layer  in  this  region.  The  heat  flux  occurs  pri¬ 
marily  in  winter  and  analysis  shows  that  much  of  the  enhancement  is  due  to 
deep  convection  which  brings  up  warm  water  and  prevents  ice  formation  in 
early  winter.  This  convection  explains  the  absence  of  ice  in  the  full 
model  simulation  far  from  the  ice  margin. 

Near  the  ice  margin  a  similar  physics  applies,  but  the  precise  loca¬ 
tion  of  the  ice  edge  becomes  more  sensitive  to  the  surface  salt  balance  and 
to  lateral  effects  in  the  oceanic  circulation.  The  sensitivity  simulations 
shown  in  Figure  24b,  for  example,  indicate  that  the  freshwater  flux  from 
melting  at  the  advancing  ice  edge  tends  to  seal  off  the  ice  margin  to  up¬ 
ward  oceanic  heat  flux  and  allows  a  farther  advance  than  would  otherwise 
occur.  A  sensitivity  simulation  without  lateral  motion  in  the  ocean  (Fig. 
24b)  also  tends  to  produce  a  more  excessive  edge  at  the  end  of  one  year, 
since  the  lateral  transport  of  heat  to  both  the  deep  and  upper  layers  is 
missing.  This  latter  result  demonstrates  the  difficulty  of  simulating  the 
seasonal  cycle  of  sea  ice  using  only  a  one-dimensional  mixed  layer  model. 

An  additional  interesting  feature  of  this  coupled  ice-ocean  model  is 
£he  modification  by  the  IcA  interaction  of  the  wind  stress  transferred  into 
the  ocean.  This  effect  is  particularly  pronounced  on  long  time  scales 
since  large  fluctuations  in  wind  stress  tend  to  average  out.  Analysis  of 
the  nine-month  average  wind  and  ice  stresses  at  buoy  locations  in  the  Arc¬ 
tic  Basin,  for  example,  shows  forces  due  to  the  ice  interaction  to  typical¬ 
ly  be  in  opposition  to  the  wind  stress  and  to  have  similar  magnitudes. 
However,  in  some  cases  the  ice  stress  can  combine  with  the  wind  stress  or 
be  at  right  angles  to  the  ice  drift  due  to  the  configuration  of  the  land 
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a.  February  50%  concentration  lim¬ 
its  from  the  full  coupled  ice/ocean 
model  and  for  the  ice-only  model, 
which  includes  a  fixed-depth  mixed 
layer.  The  observed  limit  is  taken 
from  Fleet  Weather  (Suitland,  Md. , 
U.S.  Navy)  charts  for  the  end  of 
February. 


b.  Simulated  0.05  m  thickness  con¬ 
tours  for  day  360  for  the  full  mod¬ 
el,  a  “motionless"  ocean  model,  and 
a  coupled  ice/ocean  model  with  no 
salt  fluxes  due  to  ice  freezing  and 
melting.  The  “motionless"  model 
and  the  no-salt  simulations  were 
run  for  only  1  year  and  were  ini¬ 
tialized  to  ice  thickness  and  ocean 
fields  at  the  end  of  the  fourth 
year  of  the  full  diagnostic  model 
s imulation. 


Figure  24.  Comparisons  of  simulations  of  ice  concentration  and  ice 
thickness  (from  Hibler  and  Bryan  1984). 


boundaries  and  the  ice  edge.  In  addition,  near  certain  boundaries  the  ice 
interaction  was  found  to  introduce  not  only  magnitude  changes  but  also  sign 
changes  in  the  curl  of  the  stress  transmitted  into  the  ocean.  Such  effects 
should  be  particularly  relevant  to  longer  term  prognostic  ocean  simula¬ 
tions,  and  emphasize  the  need  for  better  understanding  of  the  nature  of  ice 
interaction. 


CONCLUDING  REMARKS 

Sea  ice  dynamics  deals  with  the  way  momentum  and  energy  are  distri¬ 
buted  at  the  alr/ice/ocean  interface.  As  such,  it  plays  a  pivotal  role  in 
describing  alr-Bea  interaction  in  polar  regions.  Because  these  sources  of 
energy  and  momentum  are  often  external  to  the  ice,  understanding  sea  ice 
dynamics  entails  studying  the  coupled  air-sea-ice  system.  This  review  has 
shown  that  sea  ice  dynamics  can  cause  substantial  modification  to  the  air- 


sea  momentum,  heat,  and  salt  exchanges.  Examination  of  the  consequences  o: 
these  modifications  for  the  atmosphere  and  ocean  is  just  beginning.  From 
the  work  reviewed  here,  however,  it  is  evident  that  to  understand  these 
modifications  the  nonlinear  interaction  of  sea  ice  must  be  considered. 
Further  studies  are  needed  to  improve  our  understanding  of  the  interaction 
especially  near  the  ice  margin. 
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